User Tools

Site Tools


wiki:spec_distrib

Understanding Autoecology & Synecology to model species distribution

The main aim of this page is to show how to retrieve environmental information at point level and build up a table useful to identify Autoecology and Synecology of a species.

  • Autoecology: the branch of ecology which deals with individual species and their reactions to environmental factors. sources
  • Synecology : The branch of ecology which deals with the distribution, abundance, demography, and interactions between coexisting populations.

As an example we show a typical pattern of 3 plant species.



Autoecology and synecology of Pinus cembra, according to bioclimatic variables and to the distribution of Larix decidua and Picea abies (S. Casalegno et. al. 2010). ( script)

Autoecology of Solitary Tinamou

The Solitary Tinamou is a bird that is living in South America.
Points data can be download from ebird and gbif. Moreover, in Map Of Life there are dataset of presence and absence in National Parks. We will use these points to extract environment variable values as in this picture

gis.jpg

Prepare the data before import in R

cd ~/ost4sem/exercise/
wget https://www.dropbox.com/s/cner2yzve14adrg/Solitary_Tinamou.zip
unzip Solitary_Tinamou.zip

The following code focus in data download and data preparation. You may skip this part because the data has been already downloaded and unzipped.

# see the text data file 
# cd /home/user/ost4sem/exercise/Solitary_Tinamou/pointdata/
# more ebd_soltin1_relAug-2013.txt
# select the the column 21 LATITUDE and 22 LONGITUDE 
# awk -F  "\t" '{  print $21 ,  $22 }'   ebd_soltin1_relAug-2013.txt > lat_long_ebd.txt
 
# download WorldCom data 
# the WorldCom data have been already downloaded and stored in this directory /home/user/ost4sem/exercise/Solitary_Tinamou/worldclim
 
# wget http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/tmean_34_tif.zip 
# wget http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/tmax_34_tif.zip 
# wget http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/tmin_34_tif.zip 
# wget http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/alt_34_tif.zip 
# wget http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/tiles/cur/prec_34_tif.zip 
 
# unzipped the dataset
# unzip alt_34_tif.zip
# unzip prec_34_tif.zip
# unzip tmax_34_tif.zip
# unzip tmean_34_tif.zip
# unzip tmin_34_tif.zip
 
# consensus dataset http://www.earthenv.org/landcover 
# worldclim http://worldclim.org/version2 
# national park from here https://protectedplanet.net/

Before to start R navigate in the /home/user/ost4sem/exercise/Solitary_Tinamou/* directories and open the file with the different software that we present so far: openev , qgis , more for txt file. Try to use as much as possible the terminal to navigate in the directories and for opening the software

Install library R

Check if you have the following library are installed and if not installed them by starting R with sudo right in the terminal.

sudo R
install.packages("dismo")
install.packages("XML")
install.packages("maptools")
install.packages("rgeos")
install.packages("rasterVis")

Start R

From this point you can open in Rstudio following script /home/user/ost4sem/exercise/Solitary_Tinamou/script/Tinamou.R

library(raster)
library (dismo)
library (XML)
library (rgeos)
library (maptools)
library(sp)
library(rasterVis)
library(rgdal)
library(car)
 
## download point data of occurrences from the Global Biodiversity Information Facility (GBIF) dataset 
gbif_points = gbif('Tinamus' , 'solitarius' , download=T , geo=T)
gbif_points=gbif_points[!is.na(gbif_points$lat),]
 
## import the ebird points
ebird = read.table("/home/user/ost4sem/exercise/Solitary_Tinamou/pointdata/lat_long_ebd.txt" ,header = TRUE  )
 
## import a presence-absence dataset from parks
parks = read.table ("/home/user/ost4sem/exercise/Solitary_Tinamou/pointdata/Tinamus_solitarius_rowids.asc" , header = TRUE  )
 
## import expert range
tin_range=readOGR("/home/user/ost4sem/exercise/Solitary_Tinamou/shp/iucn_birds_proj.shp","iucn_birds_proj")
tin_range=spTransform(tin_range,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
 
## Make some pseudo-absences
pa=as.data.frame(coordinates(spsample(tin_range,70,type="random")))
 
## Build a combined dataset (with source, presence/absence, and coordinates)
points=rbind.data.frame(
  data.frame(src="gbif",obs=1,lat=gbif_points$lat,lon=gbif_points$lon),
  data.frame(src="ebird",obs=1,lat=ebird$LATITUDE,lon=ebird$LONGITUDE),
  data.frame(src="parks",obs=parks$presence,lat=parks$st_y,lon=parks$st_x),  
  data.frame(src="psuedo",obs=0,lat=pa$y,lon=pa$x)  
  )
## turn it into a spatial dataframe and define projection
coordinates(points)=c('lon','lat')
projection(points)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
 
## Create a combined src_presence field for easy plotting
points$type=paste(points$src,points$obs,sep="_")
 
##import a world country boundary to ground the map
World  = readShapePoly ("/home/user/ost4sem/exercise/Solitary_Tinamou/shp/world_country_admin_boundary_shapefile_with_fips_codes.shp")
projection(World)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
 
## check out the data
spplot(points,zcol="type",pch=1:5,col.regions=c("red","red","black","red","black"))+
  layer(sp.polygons( World))# ,xlim=c(-80 , -15) , ylim=c(-40,-5) ) ;
 
### uh oh, there are some records out in the ocean (typical GBIF!)
### There are several ways to handle this, one is to clip by the 'expert' range of the species.
 
## plot it again with the range
spplot(points,zcol="type",pch=1:5,col.regions=c("red","red","black","red","black"))+
  layer(sp.polygons( World))+
  layer(sp.polygons(tin_range,fill="grey"),under=T)
 
## So there are a few points just outside the range, but those in the ocean are most likely wrong.  
## Let's add the distance to the range polygon as a way to filter the observations:
## First we need a equidistant projection to do the calculation
dproj=CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs") 
points$dist=gDistance(spTransform(points,dproj),spTransform(tin_range,dproj),byid=T)[1,]
## that adds 'distance' (in meters) from each point to the polygon, check out the values:
hist(points$dist/1000,xlab="km to expert range")
## so some points are > 2000km from the range, let's drop any more than 10km
points=points[points$dist<10000,]
 
## plot it again with the range
spplot(points,zcol="type",pch=1:5,col.regions=c("red","red","black","red","black"))+
  layer(sp.polygons( World))+
  layer(sp.polygons(tin_range,fill="grey"),under=T)
## that looks better.
 
 
 
## plot it again adding psuedo-absences
spplot(points,zcol="type",pch=1:5,col.regions=c("red","red","black","red","black"))+
  layer(sp.polygons( World))+
  layer(sp.polygons(tin_range,fill="grey"),under=T)+
  layer(sp.points(pa,col="purple"))
 
## that looks better.
 
 
################################################################################
## Import some evironmental data (Climate & NPP) and align it to a common grid
env=stack(list.files(path = "/home/user/ost4sem/exercise/Solitary_Tinamou/worldclim/", pattern="*.tif$" , full.names = TRUE ))
## do some renaming for convenience
names(env)=sub("_34","",names(env))
names(env)[names(env)=="MOD17A3_Science_NPP_mean_00_12_clip"]="npp"
## set missing value in npp
NAvalue(env[["npp"]])=65535
 
## get total % forest
forest=sum(env[[grep("consensus",names(env))]])
names(forest)="forest"
## add forest into the env stack
env=stack(env,forest)
 
##  List all available environmental data
names(env)
 
## See how the points compare to altitude
var="prec12"
levelplot(env[[var]],col.regions=terrain.colors(20),margin=F)+
  layer(sp.polygons(tin_range))+
  layer(sp.points(points[points$obs==0,],col="black"))+ #add absences
  layer(sp.points(points[points$obs==1,],col="red"))    #add presences
 
 
## variable selection is tricky business and we're not going to dwell on it here...  
## we'll use the following variables
vars=c("tmean12","prec12","npp","alt","forest")
## look at the environmental variability
plot(env[[vars]])
 
## add the environmental data to each point
pointsd=extract(env[[vars]],points,sp=T)
 
## now check out some bivariate plots
## feel free to explore other variables
 
xyplot(npp~alt,groups=type,data=pointsd@data,auto.key=T)+
  layer(panel.ellipse(x,y,groups=pointsd$type,subscripts=T,level=.68))
 
xyplot(tmean12~prec12,groups=type,data=pointsd@data,auto.key=T)+
  layer(panel.ellipse(x,y,groups=pointsd$type,subscripts=T,level=.68))
 
xyplot(npp~forest,groups=type,data=pointsd@data,auto.key=T)+
  layer(panel.ellipse(x,y,groups=pointsd$type,subscripts=T,level=.68))
 
## Fit a very simple GLM to the data
m1=glm(obs~tmean12+prec12+npp+alt+forest,data=pointsd@data)
## Show the model summary
summary(m1)
 
## Use the model to predict p(presence) for all sites
pred=predict(env[[vars]],m1)
 
## Plot the predictions
levelplot(pred,col.regions=rainbow(20,start=.2,end=.9),margin=F)+
  layer(sp.polygons(tin_range,lwd=2))+
  layer(sp.points(points[points$obs==0,],col="black",cex=2,lwd=4))+ #add absences
  layer(sp.points(points[points$obs==1,],col="red",cex=2,lwd=4))    #add presences
 
 
## perhaps we should mask the predictions to the expert range
pred2=mask(pred,tin_range)
 
levelplot(pred2,col.regions=rainbow(20,start=.2,end=.9),margin=F)+
  layer(sp.polygons( World))+
  layer(sp.polygons(tin_range,lwd=2))+
  layer(sp.points(points[points$obs==0,],col="black",cex=2,lwd=4))+ #add absences
  layer(sp.points(points[points$obs==1,],col="red",cex=2,lwd=4))    #add presences
 
## quick accuracy assessment
predpoints=extract(pred,points,sp=T)
boxplot(layer~obs,data=predpoints@data,ylab="Estimated Propability of Occurrence",xlab="Observed Presence/Absence")
 
## to exclude the 'psuedo-absence' points, add this: [pointsd$type!="psuedo",]
## after pointsd@data
 
## what if we only used percent forest?
boxplot(forest~obs,data=pointsd@data,ylab="Percent Forest",xlab="Observed Presence/Absence")
 
##########################
## the polygons of the parks (and the associate presence/absence status) is available in 
## /home/user/ost4sem/exercise/Solitary_Tinamou/shp/protected_areas
## How would you use these polygons rather than the point data?
## check out the help for the extract function.

wiki/spec_distrib.txt · Last modified: 2018/01/24 16:46 (external edit)