Landscape – how to deal with it in science


Introduction
Some definitions
Metrics…
…and their calculation
Conclusion
Further Reading
References

…and their calculation

For a Meta-analysis it makes sense to store data in a table which contains the effect size(s) and all variables which will be used as moderator variables. Each study or, in most cases, each effect size is represented in one row. When trying to answer an ecological question we can often extract more than one effect size per study. When doing this, we have to pay attention when choosing the model we will use within the meta-analysis. But I will write more about this in the next blog-entry dedicated to the meta-analysis, where I will be writing about the details of the calculations to carry out this work.

The Code to prepare the data so that they can comfortably be used for this meta-analysis is separated into three steps. The first step is to calculate landscape metrics, the second step is to calculate correlations between abundance and landscape metrics and the third step is to paste together all effect sizes over several studies into one table. Each step is comprised in one of the functions ma-gis(), ma_prep() and ma_make.df() respectively. However, we will start with preparing GNU R that everything runs smoothly. A very first step is to save all abundance tables per the groups we want to work with (it is advisable to use one study as one of those groups). Than we save all of the below pasted scripts to the same directory from which we want to work. In this very same directory we need to create folders for each study with the same name that we will later use to define our functions (e.g. “study1”). With

setwd("~/path/to/working-directory")

we tell R where the functions will be found and where we saved the abundance tables.

source(~/path/to/script.R)

imports the functions included in this script. Now we can execute the steps 1 through 3 by running the respective functions. Attentions needs to paid to the fact that in ma_make.df() all desired studies must be included when specifying the function (i.e. ma_make.df(studies = c(“study1”, “study2”, …), buffer = 1000)). Also we can decide to calculate mean values over several years for which small mammals might have been collected. This defaults to

ma_gis(study = "study1", years = 'sampling years of study 1',  buffer = 'any desired buffer', thresh = 'a threshold for the %-cover that should be seen as forest', plot = 'F, 1 , 2 or 3, defaults to FALSE')
ma_prep(study = "study1", years = 'same years-variable as in ma_gis()', buffer = 'the buffer chosen in ma_gis()', idoi = 'vector of those clc-groups that should be treated as forest, defaults to c(23, 24, 25)')
ma_make.df(studies = c("names", "of", "files"), buffer = see ma_prep, means= 'logical, defaults to TRUE')

This will give us two files containing the raw landscape metrics and the correlations between species and landscape metrics per study, buffer and year of data-collection. ma_make.df() will use the correlations-files to coerce them to one big data-frame. The script depends on basic packages of Version 3.1.1 but for spatial calculations the packages sp, raster, maptools, rgeos, SDMTools, gfcanalysis, rasterVis and devtools; are used.

ma_gis <- function(study, years, buffer, thresh, plot=F){
&nbsp; # study = name of the study defined by you when grouping GPS-Coordinates.
&nbsp; # years = sampling years of the according study. Give as vector of form 'c(year.sampling1, year.sampling2, ...).
&nbsp; # buffer = the spatial range in which values should be taken for calculation in meters (e.g. 1000 or 1784).
&nbsp; # thresh = the threshold for "forest". Must be between 0 to 100 as "percent cover" that accounts for "forest".
&nbsp; # plot = 1 for vertical plotting, 2 for horizontal plotting and 3 for plotting of all used rasters in a 3*2 grid (does not work yet)
&nbsp;
# ---> function for checking whether all required packages are installed and loaded or install and load them if they are not installed or loaded.
&nbsp; pkgTest <- function(x){
&nbsp;&nbsp;&nbsp; if (!require(x, character.only = TRUE)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; install.packages(x, dep=TRUE)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(!require(x, character.only = TRUE)) stop("Package not found")
&nbsp;&nbsp;&nbsp; }
&nbsp; }
&nbsp;
&nbsp; pkgTest("sp"); pkgTest("raster"); pkgTest("maptools"); pkgTest("rgeos"); pkgTest("SDMTools"); pkgTest("gfcanalysis"); pkgTest("rasterVis"); pkgTest("devtools");
&nbsp; # install_github("dbahrdt/SDMTools@ignore_invalid") # this is a correction of the SDMTools function ClassStat() that a friend of me wrote for the php-code that is called by ClassStat(). It calculates edges only between those pixels that are part of the plot area and ommits the plot borders as edge, this does recently not happen (15.07.2014) when executing the function. Recently the border of the plotted area is also part of the edge and thus edge-density.

# ---> first time sample
&nbsp; t1 <- Sys.time()
&nbsp;
# ---> defining some conditions to run the function. &nbsp;
&nbsp; years.orig <- years
&nbsp; if(any(years<2000)){
&nbsp;&nbsp;&nbsp; n <- 2000-min(years)
&nbsp;&nbsp;&nbsp; years[years<2000] <- 2000
&nbsp;&nbsp;&nbsp; warning(paste("the smallest element of 'years' has a distance of ", n, " to the minimum GFC-year; interpretation will get less accurate, the bigger this distance \n any value for 'year' smaller than 2000 was set to 2000 and spatial metrics from GFC-data for those years are the 2000 values", sep=""), immediate.=T)
&nbsp; }
&nbsp; if(any(years>2012)){
&nbsp;&nbsp;&nbsp; m <- max(years)-2012
&nbsp;&nbsp;&nbsp; years[years>2012] <- 2012
&nbsp;&nbsp;&nbsp; warning(paste("the biggest element of 'years' has a distance of ", m, " to the maximum GFC-year; interpretation will get less accurate, the bigger this distance \n any value for 'year' bigger than 2012 was set to 2012 and spatial metrics from GFC-data for those years are the 2012 values", sep=""), immediate.=T)
&nbsp; }
&nbsp; if(is.na(study)){
&nbsp;&nbsp;&nbsp; stop("please specify a value for the study name 'study'")
&nbsp; }
&nbsp; if(is.na(years)[1]){
&nbsp;&nbsp;&nbsp; stop("please specify a value for the research years 'c(years)'")
&nbsp; }
&nbsp; if(is.na(buffer)){
&nbsp;&nbsp;&nbsp; stop("please specify a value for the area of interest 'buffer'")
&nbsp; }
&nbsp; if(is.na(thresh)){
&nbsp;&nbsp;&nbsp; stop("please specify a value for the threshold of forest cover 'thresh'")
&nbsp; }
&nbsp;
# ---> read in CORINE data; borders between the CORINE datasets are chosen as the median between two datasets, since it is more likely that the state of landscape was transformed to the state of the second year after half the years have passed, than the old state still exists.
&nbsp; if(median(years) < 1985) {
&nbsp;&nbsp;&nbsp; stop("only studies after 1985 can be considered")
&nbsp; }
&nbsp; print("read in target CORINE")
&nbsp; if(median(years) <= 1995){
&nbsp;&nbsp;&nbsp; clc <- raster("./CORINE land cover/clc90/g100_90.tif")
&nbsp; }
&nbsp; if(median(years) > 1995 & median(years) <= 2003){
&nbsp;&nbsp;&nbsp; clc <- raster("./CORINE land cover/clc00/g100_00.tif")
&nbsp; }
&nbsp; if(median(years) > 2003){
&nbsp;&nbsp;&nbsp; clc <- raster("./CORINE land cover/clc06/g100_06.tif")
&nbsp; }
&nbsp; breaks <- c(1:45)
&nbsp; cols <- clc@legend@colortable[1:45]

# ---> read in points. projection should be set to LAEA already before operating here. Grouping is carried out for convenience and to make automated calculation possible for points which are spread over several tiles of the GFC dataset
# --->
&nbsp; points <- readShapePoints(paste("./data/GPS/", study, "_gps.shp", sep=""))
&nbsp; crs(points) <- c("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
&nbsp; d <- dist(points@coords)
&nbsp; h <- hclust(d, method="ward.D2")
&nbsp; groups <- cutree(h, h=100000)
&nbsp; q <- 0; df.gfc2 <- NULL; df.clc2 <- NULL
&nbsp; for(l in 1:length(levels(as.factor(groups)))){
&nbsp;&nbsp;&nbsp; study.original <- study
&nbsp;&nbsp;&nbsp; points.original <- points
&nbsp;&nbsp;&nbsp; if(length(levels(as.factor(groups))) != 1) {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; print(paste("dataset will be split up into ", length(levels(as.factor(groups)))," smaller subsets", sep=""))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; study <- paste("sub", l, "_", study, sep="")
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; points <- data.frame(points)[as.numeric(names(groups[as.factor(groups)==l]))+1,]
&nbsp;&nbsp;&nbsp; coordinates(points) <- ~coords.x1+coords.x2
&nbsp;&nbsp;&nbsp; crs(points) <- c("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
&nbsp;&nbsp;&nbsp; ext <- c(points@bbox[[1]]-(buffer+1000), points@bbox[[3]]+(buffer+1000), points@bbox[[2]]-(buffer+1000), points@bbox[[4]]+(buffer+1000))

# ---> read in GFC data.
&nbsp;&nbsp;&nbsp; tpoly <- as(extent(ext), "SpatialPolygons")
&nbsp;&nbsp;&nbsp; crs(tpoly) <- crs(points); print("extracting target GFC from local files")
&nbsp;&nbsp;&nbsp; tgfc.e <- extract_gfc(tpoly, "./Global_Forest_Change"); print(paste("reprojecting target GFC to '", proj4string(points), "'", sep=""))
&nbsp;&nbsp;&nbsp; tgfc.p <- projectRaster(tgfc.e, crs=proj4string(points), method="ngb"); print("cropping target GFC to study area")
&nbsp;&nbsp;&nbsp; tgfc <- crop(tgfc.p, ext)
&nbsp;&nbsp;&nbsp; print(paste("printing object 'ra.", study, "_", buffer, "_gfc'", sep="")); assign(paste("ra.", study, "_", buffer, "_gfc", sep=""), tgfc, envir=.GlobalEnv)

# ---> crop clc to same extent.
&nbsp;&nbsp;&nbsp; print("cropping target CORINE to study area")
&nbsp;&nbsp;&nbsp; tclc <- crop(clc, ext)
&nbsp;&nbsp;&nbsp; print(paste("printing object 'ra.", study, "_", buffer, "_clc'", sep="")); assign(paste("ra.", study, "_", buffer, "_clc", sep=""), tclc, envir=.GlobalEnv)

# ---> crown-cover from GFC data.
&nbsp;&nbsp;&nbsp; print("calculating landscape metrics from GFC")
&nbsp;&nbsp;&nbsp; tgfc.e2 <- extract(tgfc[[1]], points, buffer=buffer)
&nbsp;&nbsp;&nbsp; weight <- NULL; lvl <- NULL; cover <- NULL; fc <- NULL; cc <- NULL
&nbsp;&nbsp;&nbsp; for(j in 1:length(tgfc.e2)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; x <- tgfc.e2[[j]]
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; lvl <- as.integer(levels(as.factor(x)))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cover <- c(cover, lvl); cover[cover<=thresh] <- NA
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for(i in lvl){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; weight <- c(weight, length(tgfc.e2[[j]][tgfc.e2[[j]]==i]))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; cc <- c(cc, rep(weighted.mean(cover, weight, na.rm=T), length(years)))
&nbsp;&nbsp;&nbsp; }

# ---> forest-cover is where tc00 has any value greater than defined in the variable 'thresh', forest-cover is calculated for each element in 'years'.
# ---> edge-density, et cetera from GFC data can also be extracted in the same sub-function.
&nbsp;&nbsp;&nbsp; years <- years-2000
&nbsp;&nbsp;&nbsp; fc00 <- tgfc[[1]]>thresh
&nbsp;&nbsp;&nbsp; gain <- crop(tgfc[[3]], ext)
&nbsp;&nbsp;&nbsp; loss <- crop(tgfc[[4]], ext)
&nbsp;&nbsp;&nbsp; id <- NULL; k <- NULL; yrs <- NULL; ed <- NULL; np <- NULL; pd <- NULL
&nbsp;&nbsp;&nbsp; for(j in 1:length(points)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; points.b <- gBuffer(points[j,], byid=FALSE, width=buffer, quadsegs=90)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; fc00.m <- mask(crop(fc00, extent(points.b)), rasterize(points.b, crop(fc00, extent(points.b)), silent=T)) # silent=T does not seem to work here?!
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; gain.m <- mask(crop(gain, extent(points.b)), rasterize(points.b, crop(gain, extent(points.b)), silent=T))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for(i in 1:length(years)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; loss.m <- mask(crop(loss<years[i] & loss>0, extent(points.b)), rasterize(points.b, crop(loss<years[i] & loss>0, extent(points.b))))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- (fc00.m-loss.m+gain.m)>0
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; id <- c(id, j+q)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; yrs <- c(yrs, years.orig[i])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; stats <- ClassStat(temp, cellsize=res(temp)[1])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; k <- c(k, stats$total.area[2]/(pi*buffer**2)*100)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ed <- c(ed, stats$edge.density[2]*10000) # edge.density is put out as m/m² by ClassStat
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; np <- c(np, stats$n.patches[2])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; pd <- c(pd, stats$patch.density[2]*10000)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(plot){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; print(paste("plotting ID ", id[2*j-1], " in ", years[i]+2000, sep=""))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; plot(temp)&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #par(mfrow=c(2, 2)); plot(fc00.m); plot(loss.m); plot(gain.m); plot(temp); par(mfrow=c(1, 1))
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; years <- years+2000
&nbsp;&nbsp;&nbsp; df.gfc <- data.frame(id=id, year=yrs, sub.grp = rep(l, length(groups[groups==l])*length(years)), crown.cover=cc, forest.cover=k, edge.density=ed, n.patches=np, patch.density=pd)
&nbsp;&nbsp;&nbsp; df.gfc2 <- rbind(df.gfc2, df.gfc)
&nbsp;&nbsp;&nbsp; print(paste("printing object '", study, "_", buffer, "_gfc_data'", sep=""));&nbsp; assign(paste(study, "_", buffer, "_gfc_data", sep=""), df.gfc, envir=.GlobalEnv)

# ---> landuse from CORINE data
# ---> edge density from CORINE data
&nbsp;&nbsp;&nbsp; print("calculating landscape metrics from CORINE")
&nbsp;&nbsp;&nbsp; landuse <- extract(tclc, points, buffer=buffer, weights=TRUE)
&nbsp;&nbsp;&nbsp; id <- NULL; k <- NULL; code <- NULL; ed <- NULL; ed.f <- NULL; n.patches <- NULL; area <- NULL
&nbsp;&nbsp;&nbsp; for(i in 1:length(points)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; x <- landuse[[i]]
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; y <- as.integer(levels(as.factor(x)))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for(j in y){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; code <- c(code, j)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; id <- c(id, i+q)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; points.b <- gBuffer(points[i,], byid=FALSE, width=buffer, quadsegs=90)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; clc.m <- mask(crop(tclc, extent(points.b)), rasterize(points.b, crop(tclc, extent(points.b))))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(plot){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; print(paste("plotting ID ", i+q, sep=""))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; col.sel <- as.integer(levels(as.factor(clc.m@data@values)))+1
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; plot(clc.m, col=cols[min(col.sel):max(col.sel)])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; stats <- ClassStat(clc.m, cellsize=res(clc)[1])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; clc.m[clc.m%in%c(23, 24, 25)] <- 99
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; stats2 <- ClassStat(clc.m, cellsize=res(clc)[1])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; k <- c(k, stats$total.area/(pi*buffer**2)*100)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ed <- c(ed, stats$edge.density*10000)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(any(stats2$class==99)==FALSE){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ed.f <- c(ed.f, 0)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } else {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ed.f <- c(ed.f, stats2$edge.density[stats2$class==99]*10000)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; n.patches <- c(n.patches, stats$n.patches)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; area <- c(area, stats$total.area)
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; ed.f <- data.frame(edgedensity.allforest = ed.f)
&nbsp;&nbsp;&nbsp; df.clc <- data.frame(id=as.factor(id), code=as.factor(code), n.patches=n.patches, area=area, cover=k, edge.density=ed)
&nbsp;&nbsp;&nbsp; df.clc2 <- rbind(df.clc2, df.clc)
&nbsp;&nbsp;&nbsp; print(paste("printing object '", study, "_", buffer, "_corine_data'", sep="")); assign(paste(study, "_", buffer, "_corine_data", sep=""), df.clc, envir=.GlobalEnv); assign(paste(study, "_", buffer, "_ed.all.clc", sep=""), ed.f, envir=.GlobalEnv)

# ---> plotting-functions
&nbsp;&nbsp;&nbsp; if(plot){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; print("plotting GFC and CORINE of the research area")
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; col.sel <- as.integer(levels(as.factor(tclc@data@values)))+1
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #if(plot==3){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #&nbsp; par(mfrow=c(3, 2))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #&nbsp; plot(tgfc)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #&nbsp; plot(tclc, breaks=breaks, col=cols[min(col.sel):max(col.sel)], add=T)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #&nbsp; par(mfrow=c(1, 1))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #} else{ # this does not seem to work properly with the corine plot. Please look into fixing it yourself, if needed.
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(plot==1){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; par(mfrow=c(2, 1))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(plot==2){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; par(mfrow=c(1, 2))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; plot(tclc, breaks=breaks, col=cols[min(col.sel):max(col.sel)], sub="CORINE")
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; buffer1 <- gBuffer(points, byid=TRUE, width=buffer, quadsegs = 90); plot(buffer1, col=rgb(255, 128, 0, 200, max=255), add=T)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; points(points, pch=16, col="blue")
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; text(buffer1, labels=buffer1@data[,1], pos=2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; plot(tgfc[[1]], sub="Global Forest Change")
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; buffer1 <- gBuffer(points, byid=TRUE, width=buffer, quadsegs = 90); plot(buffer1, col=rgb(255, 128, 0, 200, max=255), add=T)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; points(points, pch=16, col="blue")
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; text(buffer1, labels=buffer1@data[,1], pos=2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; par(mfrow=c(1, 1))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; #}
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; q <- length(points)
&nbsp;&nbsp;&nbsp; study <- study.original
&nbsp;&nbsp;&nbsp; points <- points.original
&nbsp; }
&nbsp; if(length(levels(as.factor(groups))) != 1) {
&nbsp;&nbsp;&nbsp; print(paste("printing object '", study, "_", buffer, "_gfc_data'", sep="")); assign(paste(study, "_", buffer, "_gfc_data", sep=""), df.gfc2, envir=.GlobalEnv)
&nbsp;&nbsp;&nbsp; print(paste("printing object '", study, "_", buffer, "_corine_data'", sep="")); assign(paste(study, "_", buffer, "_corine_data", sep=""), df.clc2, envir=.GlobalEnv)
&nbsp; }

# ---> second time sample and difference
&nbsp; t2 <- Sys.time()
&nbsp; t<- t2-t1
&nbsp; print(paste("time for processing: ", round(as.numeric(t), 2), " ", attributes(t)$units, sep=""))
&nbsp; print("done")
}
</pre>
<pre>ma_prep <- function(study, years, buffer, idoi = c(23, 24, 25)){
&nbsp; # study = name of the study defined by you when grouping GPS-Coordinates.
&nbsp; # years = sampling years of the according study. Give as vector of form 'c(year.sampling1, year.sampling2, ...).
&nbsp; # buffer = the spatial range in which values should be taken for calculation in meters (e.g. 1000 or 1784).
&nbsp; # idoi = 'id of interest'. The CORINE ID (GRID_CODE) for which edge density should be caluclated. defaults to c(23, 24, 25) since those are "true" forests.
&nbsp;
&nbsp; #library(pascal)
&nbsp; clc <- eval(as.name(paste(study, "_", buffer, "_corine_data", sep="")))
&nbsp; gfc <- eval(as.name(paste(study, "_", buffer, "_gfc_data", sep="")))
&nbsp; ed.f <- eval(as.name(paste(study, "_", buffer, "_ed.all.clc", sep="")))

# ---> read in small mammal data values
&nbsp; for(i in 1:length(years)){
&nbsp;&nbsp;&nbsp; sm.data <- read.csv(paste("./data/", study, "/", study, "_abundance_", years[i], ".csv", sep=""), row.names = 1)
&nbsp;&nbsp;&nbsp; assign(paste("sm.data", years[i], sep=""), sm.data)
&nbsp; }

# ---> set all used variables to NULL
&nbsp; evenness <- NULL; patch.number <- NULL; patch.richness <- NULL; p_i2 <- NULL; ed.all <- NULL; names <- NULL

# ---> create matrix
&nbsp; m <- matrix(data = 0, nrow=length(levels(clc$id)), ncol=length(levels(clc$code)), dimnames=list(c(levels(clc$id)), c(sort(levels(clc$code)))))

# ---> edge density (COR)
&nbsp; for(i in 1:length(idoi)){
&nbsp;&nbsp;&nbsp; ed <- NULL
&nbsp;&nbsp;&nbsp; for(k in levels(clc$id)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if((idoi[i] %in% clc$code[clc$id==k])==FALSE){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ed <- c(ed, 0)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } else {
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ed <- c(ed, round(clc$edge.density[clc$id == k & as.numeric(levels(clc$code))[clc$code] == idoi[i]], 2))&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; names <- c(names, paste("ed.", idoi[i], sep=""))
&nbsp;&nbsp;&nbsp; ed.all <- cbind(ed.all, ed); colnames(ed.all) <- names
&nbsp; }

&nbsp; for(i in levels(clc$id)){

# ---> fill matrix with available values
&nbsp;&nbsp;&nbsp; for(j in levels(clc$code)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for(n in 1:length(clc$id)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(clc$id[n] == i && clc$code[n] == j){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; m[i, j] <- round(clc[n, 5], 2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp; &nbsp;
# ---> patch number and patch richness
&nbsp;&nbsp;&nbsp; richness1 <- sum(clc$n.patches[clc$id==i])
&nbsp;&nbsp;&nbsp; patch.number <- c(patch.number, richness1)
&nbsp;&nbsp;&nbsp; richness2 <- length(clc$code[clc$id==i])
&nbsp;&nbsp;&nbsp; patch.richness <- c(patch.richness, richness2)

# ---> pielou's evenness
&nbsp;&nbsp;&nbsp; for(k in as.numeric(levels(clc$code))[clc$code][clc$id==i]){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; p_i <- sum(clc$area[clc$id==i & clc$code==k])/sum(clc$area[clc$id==i])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; p_i2 <- c(p_i2, p_i*log(p_i))
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; evenness <- c(evenness, (sum(p_i2)*-1)/log(richness1)); evenness[evenness==Inf] <- 1
&nbsp; }

# ---> add missing land cover types for comparable matrix amongst studies
&nbsp; n <- data.frame(m)
&nbsp; colnames(n) <- as.integer(colnames(m))
&nbsp; all.names <- c(1, 2, 3,&nbsp; 4, 5,&nbsp; 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44)
&nbsp; for(i in all.names){
&nbsp;&nbsp;&nbsp; if(i %in% as.integer(colnames(m))==FALSE){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; r <- rep(0, length(m[,1]))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; n <- cbind(n, r)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; colnames(n)[length(names(n))] <- i
&nbsp;&nbsp;&nbsp; }
&nbsp; }

# ---> define groups of land cover types of interest (if you want to use other groups, modify here!)
&nbsp; decid <- n[colnames(n)==23]
&nbsp; conif <- n[colnames(n)==24]
&nbsp; mixed <- n[colnames(n)==25]
&nbsp; forest_cor <- n[colnames(n)==23]+n[colnames(n)==24]+n[colnames(n)==25]
&nbsp; artif <- n[colnames(n)==1]+n[colnames(n)==2]+n[colnames(n)==3]+n[colnames(n)==4]+n[colnames(n)==5]+n[colnames(n)==6]+n[colnames(n)==7]+n[colnames(n)==8]+n[colnames(n)==9]+n[colnames(n)==10]+n[colnames(n)==11]
&nbsp; pastures <- n[colnames(n)==18]
&nbsp; het.crop <- n[colnames(n)==12]+n[colnames(n)==13]+n[colnames(n)==14]+n[colnames(n)==19]+n[colnames(n)==20]+n[colnames(n)==21]+n[colnames(n)==22]
&nbsp; perm.crop <- n[colnames(n)==15]+n[colnames(n)==16]+n[colnames(n)==17]
&nbsp; scrub_herb <- n[colnames(n)==26]+n[colnames(n)==27]+n[colnames(n)==28]+n[colnames(n)==29]
&nbsp; open_noveg <- n[colnames(n)==26]+n[colnames(n)==31]+n[colnames(n)==32]+n[colnames(n)==33]+n[colnames(n)==34]
&nbsp; wetland <- n[colnames(n)==35]+n[colnames(n)==36]+n[colnames(n)==37]+n[colnames(n)==38]+n[colnames(n)==39]
&nbsp; water <- n[colnames(n)==40]+n[colnames(n)==41]+n[colnames(n)==42]+n[colnames(n)==43]+n[colnames(n)==44]

# ---> create dataframe with all metrics
&nbsp; for(i in 1:length(years)){
&nbsp;&nbsp;&nbsp; years.orig <- years[i]
&nbsp;&nbsp;&nbsp; if(years[i]<2000){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; years[i] <- 2000
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; l <- data.frame(decid, mixed, conif, forest_cor, artif, pastures, het.crop, perm.crop, scrub_herb, open_noveg, wetland, water, patch.number, patch.richness, evenness, round(ed.f[,1], 2))
&nbsp;&nbsp;&nbsp; l <- cbind(l, ed.all)
&nbsp;&nbsp;&nbsp; l <- cbind(l, subset(gfc, gfc$year==years[i]))
&nbsp;&nbsp;&nbsp; years[i] <- years.orig
&nbsp;&nbsp;&nbsp; names(l) <- c("deciduous", "mixed", "coniferous", "forest_cor", "artificial", "pastures", "het.crop", "perm.crop", "scrub_herb", "no.veg", "wetland", "water", "patch.number", "patch.richness", "evenness", "ed.all", as.character(names), "id", "year", "crown.cover", "forest.cover", "edge.density", "n.patches", "patch.density")
&nbsp;&nbsp;&nbsp; write.table(l, paste("./data/", study, "/", study, "_vars_", buffer, "_", years[i], ".csv", sep=""), sep="\t", row.names=F)

# ---> create dataframe with correlations between abundance and variables
&nbsp;&nbsp;&nbsp; abundance <- eval(as.name(paste("sm.data", years[i], sep="")))[,5:dim(sm.data)[2]]
&nbsp;&nbsp;&nbsp; o <- matrix(data=0, nrow=length(colnames(abundance)), ncol=length(colnames(l)), dimnames=list(c(colnames(abundance)), c(colnames(l))))
&nbsp;&nbsp;&nbsp; for(k in colnames(abundance)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for(j in colnames(l)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; o[k, j] <- cor(abundance[k], l[j])
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; if(length(levels(sm.data$traplayout)) > 1){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; layout <- "several"
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; warning("several traplayouts are included in this dataset. Please split the dataset into its subsets and rerun the analysis")
&nbsp;&nbsp;&nbsp; } else{
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; layout <- levels(sm.data$traplayout)
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; if(length(levels(as.factor(sm.data$trapdensity))) > 1){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; d.trap <- "several"
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; warning("several trapdensities are included in this dataset. Please split the dataset into its subsets and rerun the analysis")
&nbsp;&nbsp;&nbsp; } else{
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; d.trap <- levels(as.factor(sm.data$trapdensity))
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; o <- data.frame(o, n_sites = cbind(length(rownames(l))), species = colnames(abundance), n.trap = sum(sm.data$trapnumber), trapnights = sum(sm.data$trapnights), layout = layout, d.trap = d.trap)
&nbsp;&nbsp;&nbsp; write.table(o, paste("./data/", study, "/", study, "_cors_", buffer, "_", years[i], ".csv", sep=""), sep="\t", row.names=F)
&nbsp; }
}
ma_make.df <- function(studies, buffer, means=T){
&nbsp; library(abind)
&nbsp; df <- NULL
&nbsp; for(i in 1:length(studies)){
&nbsp;&nbsp;&nbsp; temp <- NULL; temp2 <- NULL; temp3 <- NULL; sampleyear <- NULL; rn <- NULL; sm.data2 <- NULL
&nbsp;&nbsp;&nbsp; dir <- paste(getwd(), "/data/", studies[i], "/", sep="")
&nbsp;&nbsp;&nbsp; subsets <- dir(dir, pattern=paste(studies[i], "_cors_", buffer, sep=""))
&nbsp;&nbsp;&nbsp; s <- strsplit(studies[i], " ")[[1]]
&nbsp;&nbsp;&nbsp; name <- paste(toupper(substring(s, 1,1)), substring(s, 2), sep="", collapse=" ")
&nbsp;&nbsp;&nbsp; for(j in 1:length(subsets)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sm.data <- read.table(paste("./data/", studies[i], "/", subsets[j], sep=""), header=T)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; x <- length(colnames(sm.data))-6
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; y <- length(colnames(sm.data))-5
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; z <- length(colnames(sm.data))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sm.data2 <- sm.data[,y:z]
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(means){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(is.null(temp)){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- sm.data[,1:x]
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp[is.na(temp)] <- 0
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- data.matrix(temp)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- as.vector(temp)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp2 <- sm.data[,1:x]
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rn <- rownames(temp2); cn <- colnames(temp2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp2[is.na(temp2)] <- 0
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp2 <- data.matrix(temp2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp2 <- as.vector(temp2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- rowMeans(cbind(temp, temp2)) # does anybody know a better way to calculate means over two dataframes?
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sampleyear <- c(sampleyear, substring(subsets[j], nchar(subsets[j])-7, nchar(subsets[j])-4))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; if(length(sampleyear) > 1){
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sampleyear <- paste(min(sampleyear), " - ", max(sampleyear), sep="")
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- as.data.frame(matrix(temp, nrow=length(rn), dimnames=list(rn, cn)))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; } else{
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp2 <- sm.data[,1:x]
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp2[is.na(temp2)] <- 0
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; temp <- rbind(temp, temp2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; rn <- c(rn, rownames(temp2)); cn <- colnames(temp2)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; sampleyear <- c(sampleyear, rep(substring(subsets[j], nchar(subsets[j])-7, nchar(subsets[j])-4), length(rownames(temp2))))
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; }
&nbsp;&nbsp;&nbsp; temp3 <- cbind(temp, sm.data2, sampleyear = sampleyear, study = substr(name, 1, nchar(studies[i])-4), year = substr(studies[i], nchar(studies[i])-3, nchar(studies[i])))
&nbsp;&nbsp;&nbsp; df <- rbind(df, temp3)
&nbsp; }
&nbsp; write.table(df, paste("./data/meta-analysis", "_", buffer, ".csv", sep=""), sep="\t", row.names=FALSE)
}

The details on how the meta-analysis is calculated will follow in the next blog-post.

top

next page >>

Advertisements

One thought on “Landscape – how to deal with it in science

  1. Pingback: Landscape – how to deal with it in science | steffen ehrmann, ecologist

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s