#### Bearded Vulture Spatial analysis #### #### Aim to identify which land class (Habitat) is selected and avoided by Bearded Vulture #### ### The following steps are usually taken: ### (1) Occurrence data of the species will be compiled; ### (2) Pseudo-absence point will be created (3x the number of fixes) ### (3) Values of environmental predictor variables (i.e. Topo, Land cover) at these locations Will be extracted for each fixes ### (4) Tabulate all the results ### (5) Carry out a Data analysis (GLMM) #### Uploading and transforming Bearded Vulture restaurant data #### #### Package required #### install.packages(c('raster', 'rgdal', 'sp', 'adehabitatHR', 'rgeos', 'spatstat', 'lattice', 'FNN', 'spatialEco')) library(adehabitatHR) # package needed for home range estimation library(sp) # package needed to import and manipulate raster data-sets library(rgdal) # package needed to import ascii files into R library(raster) # package needed to manipulate raster files library(readxl) # package needed to read excel files library(rgeos) # For creating buffer library(spatstat) # Analyzing spatial point patterns library(lattice) # For visual graph library(spatialEco) # spatial data manipulation, query, sampling and modeling library(FNN) # Calculating minimum distance library(zoom) # To zoom on plots library(tibble) #### Uploading data (Fixes) #### beard = read_excel("Adult Bearded Vulture Data.xlsx") # Adult Bearded vulture tracking data NEST = read_excel("BV tracked birds nest location.xlsx") # Nesting location of Adult Bearded Vulture bv = beard[beard$Year %in% 2013:2017,] # bv setting data for year of interest (2013 to 2017) bv = bv[bv$Name != "Camo"&bv$Name != "Kloutjie",] # Removing recently adult vulture since they have no nest rm(beard) # Tidy summary(bv) str(bv) #### Projection #### prj = "+proj=longlat +datum=WGS84 +zone=35s" utm = "+proj=utm +datum=WGS84 +zone=35s" #### Visualizing tracking data and nest location #### ## Lefuma Lef = bv[bv$Name == "Lefuma",] Fnest = NEST[NEST$Bird == "Lefuma",] F13 = Lef[Lef$Year == 2013,] F14 = Lef[Lef$Year == 2014,] fx = mean(Fnest$X) # Will be used for zoom purpose fy = mean(Fnest$Y) plot(F13$X, F13$Y, col = "red", pch = 1) points(Fnest$X, Fnest$Y, pch = 3, col = "black") text(Fnest$X, Fnest$Y, label = Fnest$`SITE NAME`, cex = 1) zoomplot.zoom(fact = 5, x = fx, y = fy) plot(F14$X, F14$Y, col = "red", pch = 1) points(Fnest$X, Fnest$Y, pch = 3, col = "black") text(Fnest$X, Fnest$Y, label = Fnest$`SITE NAME`, cex = 1) zoomplot.zoom(fact = 5, x = fx, y = fy) rm(F13, F14, Fnest, Lef, fx, fy) # Tidy ## Jeremia Jer = bv[bv$Name == "Jeremia",] Jnest = NEST[NEST$Bird == "Jeremia",] jx = mean(Jnest$X) # Will be used for zoom purpose jy = mean(Jnest$Y) J13 = Jer[Jer$Year == 2013,] J14 = Jer[Jer$Year == 2014,] J15 = Jer[Jer$Year == 2015,] J16 = Jer[Jer$Year == 2016,] J17 = Jer[Jer$Year == 2017,] plot(J13$X, J13$Y, col = "red", pch = 1) points(Jnest$X, Jnest$Y, pch = 3, col = "black") text(Jnest$X, Jnest$Y, label = Jnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 5, x = jx, y = jy) plot(J14$X, J14$Y, col = "red", pch = 1) points(Jnest$X, Jnest$Y, pch = 3, col = "black") text(Jnest$X, Jnest$Y, label = Jnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 5, x = jx, y = jy) plot(J15$X, J15$Y, col = "red", pch = 1) points(Jnest$X, Jnest$Y, pch = 3, col = "black") text(Jnest$X, Jnest$Y, label = Jnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 5, x = jx, y = jy) plot(J16$X, J16$Y, col = "red", pch = 1) points(Jnest$X, Jnest$Y, pch = 3, col = "black") text(Jnest$X, Jnest$Y, label = Jnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 5, x = jx, y = jy) plot(J17$X, J17$Y, col = "red", pch = 1) points(Jnest$X, Jnest$Y, pch = 3, col = "black") text(Jnest$X, Jnest$Y, label = Jnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 5, x = jx, y = jy) rm(J13, J14, J15, J16, J17, Jer, Jnest, jx, jy) # Tidy ## Sphinx Spi = bv[bv$Name == "Sphinx",] Snest = NEST[NEST$Bird == "Sphinx",] sx = mean(Snest$X) # Will be used for zoom purpose sy = mean(Snest$Y) plot(Spi$X, Spi$Y, col = "red", pch = 1) points(Snest$X, Snest$Y, pch = 3, col = "black") zoomplot.zoom(fact = 5, x = mean(Snest$X), y = mean(Snest$Y)) text(Snest$X, Snest$Y, label = Snest$`SITE NAME`, cex = 1, pos = 4) rm(Spi, Snest, sx, sy) # Tidy ## Springbok Spr = bv[bv$Name == "Springbok",] Srnest = NEST[NEST$Bird == "Springbok",] srx = mean(Srnest$X) # Will be used for zoom purpose sry = mean(Srnest$Y) S13 = Spr[Spr$Year == 2013,] plot(S13$X, S13$Y, col = "red", pch = 1) points(Srnest$X, Srnest$Y, pch = 3, col = "black") zoomplot.zoom(fact = 5, x = srx, y = sry) text(Srnest$X, Srnest$Y, label = Srnest$`SITE NAME`, cex = 1, pos = 4) rm(Spr, Srnest, srx, sry, S13) # Tidy ## Lehlwa Leh = bv[bv$Name == "Lehlwa",] Lnest = NEST[NEST$Bird == "Lehlwa",] lx = mean(Lnest$X) # Will be used for zoom purpose ly = mean(Lnest$Y) l14 = Leh[Leh$Year == 2014,] l15 = Leh[Leh$Year == 2015,] l16 = Leh[Leh$Year == 2016,] l17 = Leh[Leh$Year == 2017,] plot(l14$X, l14$Y, col = "red", pch = 1) points(Lnest$X, Lnest$Y, pch = 3, col = "black") text(Lnest$X, Lnest$Y, label = Lnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 7, x = lx, y = ly) plot(l15$X, l15$Y, col = "red", pch = 1) points(Lnest$X, Lnest$Y, pch = 3, col = "black") text(Lnest$X, Lnest$Y, label = Lnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 7, x = lx, y = ly) plot(l16$X, l16$Y, col = "red", pch = 1) points(Lnest$X, Lnest$Y, pch = 3, col = "black") text(Lnest$X, Lnest$Y, label = Lnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 7, x = lx, y = ly) plot(l17$X, l17$Y, col = "red", pch = 1) points(Lnest$X, Lnest$Y, pch = 3, col = "black") # for 2017 nest is known look at the metadata sonja sent text(Lnest$X, Lnest$Y, label = Lnest$`SITE NAME`, cex = 1, pos = 4) zoomplot.zoom(fact = 7, x = lx, y = ly) rm(Leh, Lnest, lx, ly, l14, l15, l16, l17) # Tidy #### Creating data frame containing fixes and their respective nest locations #### ## Nest site coordinates lf = NEST[NEST$`SITE NAME` == "Witzieshoek A",] # Lefuma ln = NEST[NEST$`SITE NAME` == "Paneng B (iNkosi)",] # Inkosi Yeentaka ph = NEST[NEST$`SITE NAME` == "Column Pyramid Pass A & B (Pharoa's nest)",] # Pharoah mo = NEST[NEST$`SITE NAME` == "Kwantabanyama",] # Mollie ma = NEST[NEST$`SITE NAME` == "Khubelu",] # Mac lA = NEST[NEST$`SITE NAME` == "Sani River (Lehlwa A)",] # Lehlwa 2015 to 2017 lB = NEST[NEST$`SITE NAME` == "Sani River (Lehlwa B)",] # Lehlwa 2014 sr = NEST[NEST$`SITE NAME` == "Botelo Motlamisi C (Springbok)",] # Springbok sh = NEST[NEST$`SITE NAME` == "The Sphinx A (Sphinx)",] # Sphinx jeA = NEST[NEST$`SITE NAME` == "Amphitheatre A (Towers)",] # Jeremia Nest used for the year 2014, 2015 & 2017 jeB = NEST[NEST$`SITE NAME` == "Amphitheatre B (Ribbon Falls)",] # Jeremia Nest used for the year 2016 jeC = NEST[NEST$`SITE NAME` == "Amphitheatre C (Bilangil Falls=Jeremia's nest)",] # Jeremia Nest used for the year 2013 ## Assigning respective coordinates for each birds bv$nest_Y= ifelse(bv$Name == "Lefuma", lf$Y, ifelse(bv$Name == "Pharoah", ph$Y , ifelse(bv$Name == "Sphinx", sh$Y, ifelse(bv$Name == "Mollie", mo$Y, ifelse(bv$Name == "Mac", ma$Y, ifelse(bv$Name == "Springbok", sr$Y, ifelse(bv$Name == "InkosiYeentaka", ln$Y, ifelse(bv$Name == "Lehlwa" & bv$Year == 2015, lA$Y, ifelse(bv$Name == "Lehlwa" & bv$Year == 2016, lA$Y, ifelse(bv$Name == "Lehlwa" & bv$Year == 2017, lA$Y, ifelse(bv$Name == "Lehlwa" & bv$Year == 2014, lB$Y, ifelse(bv$Name == "Jeremia"& bv$Year == 2014, jeA$Y, ifelse(bv$Name == "Jeremia"& bv$Year == 2015, jeA$Y, ifelse(bv$Name == "Jeremia"& bv$Year == 2017, jeA$Y, ifelse(bv$Name == "Jeremia"& bv$Year == 2016, jeB$Y, ifelse(bv$Name == "Jeremia"& bv$Year == 2013, jeC$Y,"na")))))))))))))))) bv$nest_Y=as.numeric(as.character(bv$nest_Y)) bv$nest_X= ifelse(bv$Name == "Lefuma", lf$X, ifelse(bv$Name == "Pharoah", ph$X, ifelse(bv$Name == "Sphinx", sh$X, ifelse(bv$Name == "Mollie", mo$X, ifelse(bv$Name == "Mac", ma$X, ifelse(bv$Name == "Springbok", sr$X, ifelse(bv$Name == "InkosiYeentaka", ln$X, ifelse(bv$Name == "Lehlwa" & bv$Year == 2015, lA$X, ifelse(bv$Name == "Lehlwa" & bv$Year == 2016, lA$X, ifelse(bv$Name == "Lehlwa" & bv$Year == 2017, lA$X, ifelse(bv$Name == "Lehlwa" & bv$Year == 2014, lB$X, ifelse(bv$Name == "Jeremia"& bv$Year == 2014, jeA$X, ifelse(bv$Name == "Jeremia"& bv$Year == 2015, jeA$X, ifelse(bv$Name == "Jeremia"& bv$Year == 2017, jeA$X, ifelse(bv$Name == "Jeremia"& bv$Year == 2016, jeB$X, ifelse(bv$Name == "Jeremia"& bv$Year == 2013, jeC$X,"na")))))))))))))))) bv$nest_X=as.numeric(as.character(bv$nest_X)) plot(bv$nest_X, bv$nest_Y, cex=0.8, pch=16, col="purple") # Visualizing nest locations rm(jeA, jeB, jeC, lf, lA, lB, ln, ma, mo, ph, sh, sr) # Tidy #### finding fixes distance from nest for each birds #### ## Add distance to nest using cartesian equation: mdist = function(x,y,x2,y2){ sqrt((x-x2)^2+(y-y2)^2) } ## transforming fixes bv2 = bv # making a copy of lat-long data coordinates(bv2) = ~ X + Y # Project fix location crs(bv2) = CRS(prj) # Giving current projection bv2 = spTransform(bv2, CRS(utm)) # Transforming to UTM for measurement bv2 = as.data.frame(bv2) # Return to data-frame ## transforming nest location coordinates(bv2) = ~ nest_X + nest_Y # same procedure as above crs(bv2) = CRS(prj) bv2 = spTransform(bv2, CRS(utm)) bv2 = as.data.frame(bv2) head(bv2) ## apply mdist function (calculate the distances from fixes to nest) bv$nest_dist = mdist(bv2$X, bv2$Y, bv2$nest_X, bv2$nest_Y)/1000 summary(bv$nest_dist) rm(mdist, NEST) ## Finding the home range for each bird and specific years quantile(bv$nest_dist, probs = 0.99) #### remove outlier points using home range circle buffer #### bv = subset(bv, bv$nest_dist <= 95.54541) bv = bv[,-14] # remove the column because will re-calculate the nest distance together with pseudo points rm(bv2) # Tidy #### Creating random point (Pseudo-absence points) #### buff = 95.54 #km pseudo_1 = NULL ID = unique(bv$Name) for(i in 1:length(ID)){ bird = bv[bv$Name == ID[i],] # bv data for one individuals nest = bird[1:1,] # Extracting first line of data for that ID coordinates(nest)=~ nest_X + nest_Y # make nest coordinate into spatial points crs(nest)= crs(prj) # set coordinate ref system nest = spTransform(nest, crs(utm)) # transform to utm so we can measure buffer size in meter buffer = gBuffer(nest, width = buff*1000, byid=F, capStyle = "ROUND") # create buffer buffer = spTransform(buffer, crs(prj)) # transform buffer back to long-lat so the random points will be in the same coord-sys ps_points = spsample(buffer,n=(length(bird$Name)*3),type="random") #create random points. 3 times the number of real points per bird ps_points = as.data.frame(ps_points) # transform to data frame ps_points$name = c(as.character(bird$Name)) # set classes ps_points$nest_Y = bird$nest_Y[1] ps_points$nest_X = bird$nest_X[1] names(ps_points)[names(ps_points) == "y"] <- "Y" # names to match real points data-frame names(ps_points)[names(ps_points) == "x"] <- "X" plot(ps_points$X, ps_points$Y, col="red") points(bird$X, bird$Y, col="blue") pseudo_1 = rbind(pseudo_1, ps_points) } nrow(pseudo_1)/3==nrow(bv) #TRUE check there are the right number of points rm(bird, nest, ps_points, buff, i, ID) # Tidy ### Creating random point EQUAL NUMBER (Pseudo-absence points) ### buff = 95.54 #km pseudo_2 = NULL ID = unique(bv$Name) for(i in 1:length(ID)){ bird = bv[bv$Name == ID[i],] # bvset data for one individuals nest = bird[1:1,] # Extracting first line of data for that ID coordinates(nest)=~ nest_X + nest_Y # make nest coordinate into spatial points crs(nest)= crs(prj) # set coordinate ref system nest = spTransform(nest, crs(utm)) # transform to utm so we can measure buffer size in meter buffer = gBuffer(nest, width = buff*1000, byid=F, capStyle = "ROUND") # create buffer buffer = spTransform(buffer, crs(prj)) # transform buffer back to long-lat so the random points will be in the same coord-sys ps_points = spsample(buffer,n=(length(bird$Name)),type="random") #create random points EQUAL number of real points per bird ps_points = as.data.frame(ps_points) # transform to data frame ps_points$name = c(as.character(bird$Name)) # set classes ps_points$nest_Y = bird$nest_Y[1] ps_points$nest_X = bird$nest_X[1] names(ps_points)[names(ps_points) == "y"] <- "Y" # names to match real points data-frame names(ps_points)[names(ps_points) == "x"] <- "X" plot(ps_points$X, ps_points$Y, col="red") points(bird$X, bird$Y, col="blue") pseudo_2 = rbind(pseudo_2, ps_points) } nrow(pseudo_2)==nrow(bv) #TRUE check there are the right number of points rm(bird, nest, ps_points, buff, i, ID) # Tidy ## columns names need to match to be able to join data-frames: names(pseudo_1) names(pseudo_2) names(bv) ## Cleaning up bearded vulture data-frame bv$Date = NULL bv$Time = NULL bv$Speed = NULL bv$Course = NULL bv$Age = NULL ## Adding column to pseudo absence data-frame pseudo_1$Sex = "na" pseudo_1$Year = "na" pseudo_1$Altitude = "na" pseudo_2$Sex = "na" pseudo_2$Year = "na" pseudo_2$Altitude = "na" ## creating two different set of bearded vulture data-frame bv_3 = bv bv_e = bv ## add binary values to both data-frames (will be used in analyses) and join pseudo_1$bin = 0 bv_3$bin = 1 pseudo_2$bin = 0 bv_e$bin = 1 names(pseudo_1)[names(pseudo_1) == "name"] = "Name" # Output should be TURE names(pseudo_2)[names(pseudo_2) == "name"] = "Name" # Output should be TURE bv_3 = rbind(bv_3, pseudo_1) bv_e = rbind(bv_e, pseudo_2) rm(pseudo_1, pseudo_2, buffer, bv) # Tidy #### adding distance of fixes and pseudo absence point from nest #### ## Add distance to nest using cartesian equation: mdist = function(x,y,x2,y2){ sqrt((x-x2)^2+(y-y2)^2) } ## transforming fixes for data-frame 3 x pseudo bv23 = bv_3 # making a copy of lat-long data coordinates(bv23) = ~ X + Y # Project fix location crs(bv23) = CRS(prj) # Giving current projection bv23 = spTransform(bv23, CRS(utm)) # Transforming to UTM for measurement bv23 = as.data.frame(bv23) # Return to data-frame ## transforming fixes for data-frame EQUAL pseudo bv2e = bv_e # making a copy of lat-long data coordinates(bv2e) = ~ X + Y # Project fix location crs(bv2e) = CRS(prj) # Giving current projection bv2e = spTransform(bv2e, CRS(utm)) # Transforming to UTM for measurement bv2e = as.data.frame(bv2e) # Return to data-frame ## transforming nest location for data-frame 3 x pseudo coordinates(bv23) = ~ nest_X + nest_Y # same procedure as above crs(bv23) = CRS(prj) bv23 = spTransform(bv23, CRS(utm)) bv23 = as.data.frame(bv23) ## transforming nest location for data-frame EQUAL pseudo coordinates(bv2e) = ~ nest_X + nest_Y # same procedure as above crs(bv2e) = CRS(prj) bv2e = spTransform(bv2e, CRS(utm)) bv2e = as.data.frame(bv2e) head(bv23) class(bv23) ## apply minimum dist function (calculate the distances from fixes to nest) bv_3$nest_dist = mdist(bv23$X, bv23$Y, bv23$nest_X, bv23$nest_Y)/1000 bv_e$nest_dist = mdist(bv2e$X, bv2e$Y, bv2e$nest_X, bv2e$nest_Y)/1000 rm(mdist, bv23, bv2e) # Tidy #### Nearest distance to vulture restaurant #### ## bv setting data of interest (area around Lesotho) vulr = read_excel("Vulture restaurant data.xlsx") # Uploading VR dataset province = c("KwaZulu Natal", "KwaZulu-Natal", "KZN", "Eastern Cape", "Free State") # selecting Province of interest sfs = subset(vulr, vulr$Province %in% province) sfs = subset(sfs, sfs$Status_Category %in% "Active") sfs = subset(sfs, sfs$Verified_2018 %in% "Y") sfs = subset(sfs, sfs$Structure %in% "SFS") sfs = subset(sfs, select = c('SFS_Code', 'Restaurant Name', 'Y', 'X', 'Provisioning_Category')) rm(province, vulr) # Tidy ## filtering vulture restaurant falling outside buffer zones sfs2 = sfs # Making copy of original data-set sfs2$Y = as.numeric(as.character(sfs2$Y)) sfs2$X = as.numeric(as.character(sfs2$X)) coordinates(sfs2) = ~ X + Y # Project fix location crs(sfs2) = CRS(prj) # Giving current projection NEST = read_excel("BV tracked birds nest location.xlsx") crs(NEST) # Checking nest data projection coordinates(NEST) = ~ X + Y crs(NEST) = crs(prj) Nest = spTransform(NEST, crs(utm)) # Projecting in UTM to calculate buffer in meter buff = 95.5 bufpoint = gBuffer(Nest, width = buff*1000, capStyle = "ROUND", byid=F) # created buffer for all nest by id should be FALSE otherwise all the polygon get merge bufpoint = spTransform(bufpoint, crs(prj)) plot(bufpoint) vr = crop(sfs2, bufpoint) # Crop point that fall inside the buffers vr = as.data.frame(vr) # Create a data-frame vr$Provisioning_Category = as.factor(as.character(vr$Provisioning_Category)) vr$Provisioning_Category # should have six categories write.csv(vr,"vr_filtered_dataset.csv", row.names = F, col.names = T) NEST = as.data.frame(NEST) # change from spatial coordinate to data frame rm(sfs2,Nest,sfs,buff, bufpoint, NEST) # Tidy ### calculating minimum distance from fixes + pseudo points to vulture restaurant ### vr = read.csv("vr_filtered_dataset.csv") vr2 = vr coordinates(vr2) =~ X + Y crs(vr2)= crs(prj) vr2 = spTransform(vr2, crs(utm)) bv23 = bv_3 # Data-set with 3 x pseudo points coordinates(bv23) =~ X + Y crs(bv23)= crs(prj) bv23 = spTransform(bv23, crs(utm)) bv2e = bv_e # Data-set with EQUAL pseudo points coordinates(bv2e) =~ X + Y crs(bv2e)= crs(prj) bv2e = spTransform(bv2e, crs(utm)) s_dist3 = apply(gDistance(bv23, vr2,byid=T),2,min)/1000 # calculate minimum distance (/1000 to make km) s_diste = apply(gDistance(bv2e, vr2,byid=T),2,min)/1000 # calculate minimum distance (/1000 to make km) bv_3$vr_dist = s_dist3 bv_e$vr_dist = s_diste bv_3$Near_ID = as.vector(apply(gDistance(bv23, vr2,byid=T), 2, function(x) which(x==min(x)))) # assigning ID to restaurant near to the points bv_3$Near_ID = as.factor(as.character(bv_3$Near_ID)) bv_e$Near_ID = as.vector(apply(gDistance(bv2e, vr2,byid=T), 2, function(x) which(x==min(x)))) # assigning ID to restaurant near to the points bv_e$Near_ID = as.factor(as.character(bv_e$Near_ID)) ## Assigning provisioning categories to the nearest vr ## bv_3$Near_ID = as.character(bv_3$Near_ID) # Converting to character vl = c(4,5,6,7,8,10,13,15,19,21,23,25,28,31,32,33,35,38,39,40,41,42,43,49,50,51,53,55,56) l = c(1,9,11,12,18,30,37,47,48,52,54) m = c(2,14,16,17,20,22,24,29,34,36,44,45) h = c(27,46) vh = 3 uk = 26 bv_3$Near_ID[bv_3$Near_ID %in% vh] = "Very High" bv_3$Near_ID[bv_3$Near_ID %in% h] = "High" bv_3$Near_ID[bv_3$Near_ID %in% m] = "Medium" bv_3$Near_ID[bv_3$Near_ID %in% l] = "Low" bv_3$Near_ID[bv_3$Near_ID %in% vl] = "Very Low" bv_3$Near_ID[bv_3$Near_ID %in% uk] = "Unknown" bv_3$Near_ID = as.factor(bv_3$Near_ID) # Converting into categorical bv_e$Near_ID = as.character(bv_e$Near_ID) bv_e$Near_ID[bv_e$Near_ID %in% vh] = "Very High" bv_e$Near_ID[bv_e$Near_ID %in% h] = "High" bv_e$Near_ID[bv_e$Near_ID %in% m] = "Medium" bv_e$Near_ID[bv_e$Near_ID %in% l] = "Low" bv_e$Near_ID[bv_e$Near_ID %in% vl] = "Very Low" bv_e$Near_ID[bv_e$Near_ID %in% uk] = "Unknown" bv_e$Near_ID = as.factor(bv_e$Near_ID) # Converting into categorical ## renaming Near_ID column ## names(bv_3)[names(bv_3) == 'Near_ID'] = 'provisioning_cat' names(bv_e)[names(bv_e) == 'Near_ID'] = 'provisioning_cat' rm(s_dist3, s_diste, bv23, bv2e, vr2, l, vl, vr, m, h, vh, uk) # Tidy #### Extracting topography data #### ## Uploading the layer ## srtm = raster("southern_africa.tif") # Topography layer ## cropping area of interest ## ext = extent(26.5, 31.5, -32, -27) # (first row: xmin, xmax; second row: ymin, ymax) dem = crop(srtm, ext) plot(dem) rm(srtm, ext) # Tidy ## creating topographic parameters ## slope = terrain(dem, opt="slope", unit='degrees', neighbors=8) # make slope layer from DEM tri = terrain(dem, opt="TRI" ) # make terrain ruggedness index roughness = terrain(dem, opt="roughness") # make roughness layer vrm=spatialEco::vrm(dem, s=5) # make vector ruggedness measure # #create terrain stack for easy extractions ## terrain.stack = stack(list(elev=dem, slope=slope, tri=tri, roughness=roughness, vrm=vrm)) ## Extraction ## bv23 = bv_3 # making copy of lat-long data coordinates(bv23) = ~ X + Y # Project fix location crs(bv23) = CRS(prj) # Giving current projection bv2e = bv_e # making copy of lat-long data coordinates(bv2e) = ~ X + Y # Project fix location crs(bv2e) = CRS(prj) # Giving current projection vars3=extract(terrain.stack, bv23) #extract vars3=as.data.frame(vars3) varse=extract(terrain.stack, bv2e) #extract varse=as.data.frame(varse) bv_3 = cbind(bv_3, vars3) bv_e = cbind(bv_e, varse) rm(dem, roughness, slope, tri, terrain.stack, vars3, varse, bv23, bv2e, vrm) # Tidy #### Extracting land cover data #### GLC = raster("glc_100m_2015.tif") # Global land cover layer 100m resolution 2016 projection(GLC) ## cropping area of interest ## ext = extent(26.5, 31.5, -32, -27) # (first row: xmin, xmax; second row: ymin, ymax) glc = crop(GLC, ext) rm(GLC, ext) # Tidy ## Reclassifying Land cover classes ## vals = unique(values(glc)) vals recl = matrix(c(vals, c(1,2,3,4,5,6,7,8, rep(5, 4),0, rep(5,3), 7)), ncol = 2) recl # 1 = Grassland, 2 = Shrubland, 3 = Cultivated land, 4 = Urban, 5 = Forest, 6 = Wetland, 7 = Waterbodies # 8 = Bareland lc = reclassify(glc, rcl = recl) plot(lc) lc = as.factor(lc) levels(lc) ## Add names of land-cover categories to raster ## land_cover = levels(lc)[[1]] land_cover[,"landcover"] = c("NoData", "Grassland" ,"Shrubland", "Cultivated", "Built_up", "Forest", "wetland", "Waterbodies", "Bareland") levels(lc) = land_cover plot(lc) levels(lc) writeRaster(lc, filename="tidy_land_cover.tif") rm(glc, recl, vals, land_cover) # Tidy ## Extraction ## bv23 = bv_3 # making copy of lat-long data coordinates(bv23) = ~ X + Y # Project fix location crs(bv23) = CRS(prj) # Giving current projection bv2e = bv_e # making copy of lat-long data coordinates(bv2e) = ~ X + Y # Project fix location crs(bv2e) = CRS(prj) # Giving current projection cover3 = extract(lc, bv23) cover3 = as.data.frame(cover3) bv_3 = cbind(bv_3, cover3) covere = extract(lc, bv2e) covere = as.data.frame(covere) bv_e = cbind(bv_e, covere) ## renaming column ## names(bv_3)[names(bv_3) == 'cover3'] = 'land_cover' names(bv_e)[names(bv_e) == 'covere'] = 'land_cover' ## renaming number to land cover code ## bv_3$land_cover = as.character(bv_3$land_cover) # Converting to character bv_3$land_cover[bv_3$land_cover == "1"] = "Gra" bv_3$land_cover[bv_3$land_cover == "2"] = "Shr" bv_3$land_cover[bv_3$land_cover == "3"] = "Cul" bv_3$land_cover[bv_3$land_cover == "4"] = "Urb" bv_3$land_cover[bv_3$land_cover == "5"] = "For" bv_3$land_cover[bv_3$land_cover == "6"] = "Wet" bv_3$land_cover[bv_3$land_cover == "7"] = "Wat" bv_3$land_cover[bv_3$land_cover == "8"] = "Bar" bv_e$land_cover = as.character(bv_e$land_cover) # Converting to character bv_e$land_cover[bv_e$land_cover == "1"] = "Gra" bv_e$land_cover[bv_e$land_cover == "2"] = "Shr" bv_e$land_cover[bv_e$land_cover == "3"] = "Cul" bv_e$land_cover[bv_e$land_cover == "4"] = "Urb" bv_e$land_cover[bv_e$land_cover == "5"] = "For" bv_e$land_cover[bv_e$land_cover == "6"] = "Wet" bv_e$land_cover[bv_e$land_cover == "7"] = "Wat" bv_e$land_cover[bv_e$land_cover == "8"] = "Bar" bv_3$land_cover = as.factor(bv_3$land_cover) # Converting into categorical bv_e$land_cover = as.factor(bv_e$land_cover) # Converting into categorical rm(bv23, bv2e, cover3, covere, lc) # Tidy #### Creating aggregated urban cover layer #### urb1 = raster("urban.tif") # urban cover layer (processed in ArcGIS due to memory limitation in R) buff = 100 nest = read_excel("BV tracked birds nest location.xlsx") coordinates(nest)=~ X + Y # make nest coordinate into spatial points crs(nest)= crs(prj) # set coordinate ref system nest = spTransform(nest, crs(utm)) # transform to utm so we can measure buffer size in meter buffer = gBuffer(nest, width = buff*1000, byid=F, capStyle = "ROUND") # create buffer buffer = spTransform(buffer, crs(prj)) plot(buffer) urb = crop(urb1, buffer) urbm = mask(urb, buffer) plot(urbm) library(igraph) urban_clumps = clump(urbm, directions=8) #this looks for clumps in 8 directions - ie; each side and each cell on the diagonal corners. If you change it to directions=4 then the clumps are found using the cells directly to the side, not on the corners. clumpFreq = freq(urban_clumps) #this counts the number of cells in each clump clumpFreq = as.data.frame(clumpFreq) plotfreq = clumpFreq hist(plotfreq$count) plotfreq = subset(plotfreq, count<30) #limit xlim and plot again par(pty='s') hist(plotfreq$count, xlab = "Built-up pixel clump count", ylab = "Frequency", main = "") box() rm(buffer, clumpFreq, plotfreq, urb, urb1, urbm, buff, urban_clumps) ## Original urban layer of 30m will be projected on a coarser scale by mean so as to extract value## ## NOTE: 0.001 degree = 111m (convention used for the aggregation) ## two different aggregate will be used to see which one is better for the extraction ## ## Aggregate from 30x30 resolution to 300x300 (factor = 10) urb300 = aggregate(urb1, fact = 10, fun = mean) plot(urb300) writeRaster(urb300, filename="urban_mean_300.tif") ## Aggregate from 30x30 resolution to 210x210 (factor = 7) urb210 = aggregate(urb1, fact = 7, fun = mean) plot(urb210) writeRaster(urb210, filename="urban_mean_210.tif") rm(urb1, urb300, urb200) # Tidy #### Mean urban Extraction #### bv23 = bv_3 # making copy of lat-long data coordinates(bv23) = ~ X + Y # Project fix location crs(bv23) = CRS(prj) # Giving current projection bv2e = bv_e # making copy of lat-long data coordinates(bv2e) = ~ X + Y # Project fix location crs(bv2e) = CRS(prj) # Giving current projection ### Extraction 300m resolution ### urb_300 = raster("urban_mean_300.tif") urb_3 = extract(urb_300, bv23) urb_3 = as.data.frame(urb_3)*100 bv_3 = cbind(bv_3, urb_3) urb_e = extract(urb_300, bv2e) urb_e = as.data.frame(urb_e)*100 bv_e = cbind(bv_e, urb_e) ### Extraction 210m resolution ### urb_210 = raster("urban_mean_210.tif") urb_33 = extract(urb_210, bv23) urb_33 = as.data.frame(urb_33)*100 bv_3 = cbind(bv_3, urb_33) urb_ee = extract(urb_210, bv2e) urb_ee = as.data.frame(urb_ee)*100 bv_e = cbind(bv_e, urb_ee) rm(urb_300, urb_210, urb_e, urb_ee, urb_m_300, urb_l_300, urb_h_300, bv23, bv2e, urb_3, urb_33) # Tidy ## Renaming column ## names(bv_3)[names(bv_3) == 'urb_3'] = 'pr_urban_300' names(bv_e)[names(bv_e) == 'urb_e'] = 'pr_urban_300' names(bv_3)[names(bv_3) == 'urb_33'] = 'pr_urban_210' names(bv_e)[names(bv_e) == 'urb_ee'] = 'pr_urban_210' rm(urban3, urbane, bv23, bv2e, urb, pr_urban) # Tidy #### Measuring minimum distance from urban pixel to fixes/pseudo #### urb1 = raster("urban.tif") # urban cover layer (processed in ArcGIS due to memory limitation in R) levels(urb1) extent(urb1) nest = read_excel("BV tracked birds nest location.xlsx") # Read nest locations data crs(nest) # Checking nest data projection coordinates(nest) = ~ X + Y # Converting to spatial coordinates crs(nest) = crs(prj) # Assigning projection nest = spTransform(nest, crs(utm)) buff = 100 bufpoint = gBuffer(nest, width = buff*1000, capStyle = "ROUND", byid=T) bufpoint = spTransform(bufpoint, crs(prj)) ext = aggregate(bufpoint) # Aggregation was done to Crop area of interest easily plot(ext) rm(nest,bufpoint,buff, urb_level) # Tidy urb_c = crop(urb1, ext) # cropping urb_m = mask(urb_c, ext) # Masking (did it in two part it goes faster) recl = c(0 , NA, 1, 1) # 0 (non-urban) is converted to NA and 1 remains 1 (urban) recl_m = matrix(recl, ncol = 2, byrow = T) # Creat Matrix to be used for reclassification urban_m = reclassify(urb_m, rcl = recl_m) # reclassification writeRaster(urban_m, filename="urban_m.tif") # Create raster layer in folder rm(recl,recl_m,urb_c, urb_m, urb1, ext) # Tidy urban_m = raster("urban_m.tif") plot(urban_m) ## Get lesotho border countries = raster::getData("ISO3") countries[which(countries$NAME=="Lesotho"),] les = raster::getData('GADM',country='LSO',level=0) les rm(countries) rasterOptions() options(rasterTmpTime = 1) rasterOptions(timer = TRUE) rasterOptions(chunksize = 5e8) # 500MB ## Read in track data bv2 = bv # Bearded vulture data coordinates(bv2) =~ X + Y # Turn into spatial coordinates crs(bv2)= crs(prj) # Assigning projection bv2 = spTransform(bv2, crs(utm)) ## Turning urban pixel into spatial points urban_m = raster("urban_m.tif") urban_point = rasterToPoints(urban_m, spatial = FALSE) head(urban_point) urban_point = as.data.frame(urban_point) write.csv(urban_point,"urban_point.csv", quote = F, row.names = F , col.names = T) urban_point = read.csv("urban_point.csv", header = T) coordinates(urban_point) = ~ x + y crs(urban_point) = crs(prj) urban_point = spTransform(urban_point, crs(utm)) ## ANOTHER APPROCHE TO RUME SINGLE PIXEL TO LIMIT ERROR ## #try to remove single urban points : #could try drop_crumbs function in smoothr package instead here: library(igraph) #urban_m = raster("urban_m.tif") #urban_clumps <- clump(urban_m, directions=8) #find clumps #clumpFreq <- freq(urban_clumps) #clumpFreq <- as.data.frame(clumpFreq) #Coerce freq table to data.frame #excludeID <- clumpFreq$value[which(clumpFreq$count<=5)] #Put these into a vector of clump ID's to be removed #urban_clumps[urban_clumps %in% excludeID] <- NA #writeRaster(urban_clumps, filename="urban_clump(5).tif") #plot(urban_clumps) #looks similar to first one ## END ## ## calculate track point distances urb_dist = get.knnx(data = coordinates(urban_point), query = coordinates(bv2), k = 1, algorithm = "kd_tree") bv$mini_urbdist = urb_dist$nn.dist/1000 rm(bv2, urb_dist, urban_m, urban_point) # Tidy ## visualize in Google Earth by creating KML file urban_point = read.csv("urban_point.csv", header = T) coordinates(urban_point) = ~ x + y crs(urban_point) = crs("+proj=longlat +datum=WGS84") writeOGR(urban_point, dsn="urban.kml", layer= "urban", driver="KML") #### exporting final data set #### write.csv(bv_3,"bv_final_dataset_x3.csv", row.names = F, col.names = T) write.csv(bv_e,"bv_final_dataset_eq.csv", row.names = F, col.names = T) rm(bv_3, bv_e, prj, utm) # Tidy #### Data analysis #### #### Package required #### install.packages(c('lme4', 'MuMIn', 'ROCR', 'ggplot2', 'ggcorrplot', 'visreg', 'car', 'dplyr')) library(lme4) # For function Generalized Linear Mix Model (GLMM) library(MuMIn) # For Model selection library(ROCR) # To assess the Model performance library(ggplot2) # To create plots library(ggcorrplot) # Create Pearson correlation plot library(visreg) # Displaying the results of a fitted model library(car) # for vif and correlation analysis library(dplyr) # data exploration library(effects) # Visualization of model #### Upload data #### bv3 = read.csv("bv_final_dataset_x3.csv") # dataset with 3 time pseudo absence points bvE = read.csv("bv_final_dataset_eq.csv") # dataset with equal pseudo absence points #### Trimming down of the data-frame with variables of interest #### trim = c("nest_X", "nest_Y", "Altitude") bv3[trim] = NULL bvE[trim] = NULL rm(trim) # Tidy #### Check data in term of character, integer and factor #### str(bv3) # This show each column and type of data stored (int or character or numeric) cat = c("Year", "Name", "Sex", "land_cover") bv3[cat] = lapply(bv3[cat], factor) # as.factor() could also be used bvE[cat] = lapply(bvE[cat], factor) # as.factor() could also be used str(bv3) rm(cat) # Tidy #### Data exploration #### # Exploring the raw data is a crucial part of any scientific study. Zurr et al. (2009) suggest researchers follow a # 8 step data exploration protocol before proceeding with data analysis in order to avoid common statistical errors # (see Figure 1, pp.4, Zurr et al. 2009). For example, because my response variable (Y) is binomial, # I did not test for outliers in the response variable (e.g., Step 1: Outliers in X and Y, Zurr et al. 2009), # nor did I check for normality, homogeneity or zero trouble in Y (e.g., steps 2-4, Fig. 1, Zurr et al. 2009). # I also did not check for independence in the response variable (e.g., Step 8: Independence in Y, Zurr et al. 2009) # because I have a prior knowledge that when GPS locations are repeatedly sampled from individuals over short time intervals # (as is the case with this sample of GPS locations), they will exhibit some degree of spatial-temporal autocorrelation. # I was able to explore the data for (a) outliers in the explanatory variables; (b) collinearity between explanatory variables; # (c) relationships between the response and explanatory variables; and (d) interactions between explanatory variables. ### (a - i) exploring continuous explanatory variables ### par(mfrow=c(1,2)) # ELEVATION boxplot(bv3$elev, ylab="Elevation (m)", col="blue", cex.lab=0.8, pch=16) dotchart(bv3$elev, xlab="Elevation (m)", ylab="Order of Data", pch=16, cex.lab=0.8) boxplot(bvE$elev, ylab="Elevation (m)", col="blue", cex.lab=0.8, pch=16) dotchart(bvE$elev, xlab="Elevation (m)", ylab="Order of Data", pch=16, cex.lab=0.8) # SLOPE boxplot(bv3$slope, ylab="Slope (m)", col="green", cex.lab=0.8, pch=16) dotchart(bv3$slope, xlab="Slope (m)", ylab="Order of Data", pch=16, cex.lab=0.8) boxplot(bvE$slope, ylab="Slope (m)", col="green", cex.lab=0.8, pch=16) dotchart(bvE$slope, xlab="Slope (m)", ylab="Order of Data", pch=16, cex.lab=0.8) # TERRAIN RUGGEDNESS INDEX boxplot(bv3$tri, ylab="TRI", col="red", cex.lab=0.8, pch=16) dotchart(bv3$tri, xlab="TRI", ylab="Order of Data", pch=16, cex.lab=0.8) boxplot(bvE$tri, ylab="TRI", col="red", cex.lab=0.8, pch=16) dotchart(bvE$tri, xlab="TRI", ylab="Order of Data", pch=16, cex.lab=0.8) # ROUGHNESS boxplot(bv3$roughness, ylab="Roughness", col="yellow", cex.lab=0.8, pch=16) dotchart(bv3$roughness, xlab="Roughness", ylab="Order of Data", pch=16, cex.lab=0.8) boxplot(bvE$roughness, ylab="Roughness", col="yellow", cex.lab=0.8, pch=16) dotchart(bvE$roughness, xlab="Roughness", ylab="Order of Data", pch=16, cex.lab=0.8) # VRM boxplot(bv3$vrm, ylab="VRM", col="purple", cex.lab=0.8, pch=16) dotchart(bv3$vrm, xlab="VRM", ylab="Order of Data", pch=16, cex.lab=0.8) boxplot(bvE$vrm, ylab="VRM", col="purple", cex.lab=0.8, pch=16) dotchart(bvE$vrm, xlab="VRM", ylab="Order of Data", pch=16, cex.lab=0.8) # PROPORTION URBAN boxplot(bv3$pr_urban, ylab="proportion urban", col="grey", cex.lab=0.8, pch=16) dotchart(bv3$pr_urban, xlab="proportion urban", ylab="Order of Data", pch=16, cex.lab=0.8) boxplot(bvE$pr_urban, ylab="proportion urban", col="grey", cex.lab=0.8, pch=16) dotchart(bvE$pr_urban, xlab="proportion urban", ylab="Order of Data", pch=16, cex.lab=0.8) ### (a - ii) exploring categorical explanatory variables ### par(mfrow=c(1,2)) # LAND COVER barplot(table(bv3$land_cover), main="Distribution of Points Between Habitat Classes", xlab ="Habitat Class", ylab="Count") barplot(table(bvE$land_cover), main="Distribution of Points Between Habitat Classes", xlab ="Habitat Class", ylab="Count") # VULTURE RESTAURANT barplot(table(bv3$provisioning_cat), main="Distribution of Points Between VR", xlab ="provisioning category", ylab="Count") barplot(table(bvE$provisioning_cat), main="Distribution of Points Between VR", xlab ="provisioning category", ylab="Count") ### (b - i) collinearity between continuous explanatory variables ### cor3 = cor(bv3[,c(6, 7, 8, 9, 10, 11, 12, 14, 15)], method = "pearson") round(cor3, 2) ggcorrplot(cor3, hc.order = TRUE, type = "lower", lab = TRUE) corE = cor(bvE[,c(6, 7, 8, 9, 10, 11, 12, 14, 15)], method = "pearson") round(corE, 2) ggcorrplot(corE, hc.order = TRUE, type = "lower", lab = TRUE) rm(cor3, corE) # Tidy ### variance inflation factor ### vm1 = glm(bin ~ elev + slope + vrm + tri + roughness + pr_urban + nest_dist + vr_dist, data = bv3, family = binomial()) vm2 = glm(bin ~ elev + slope + vrm + tri + roughness + pr_urban + nest_dist + vr_dist, data = bvE, family = binomial()) vif(vm1) vif(vm2) # slope, tri and roughness has the highest value confirming the correlation graph vm3 = glm(bin ~ elev + vrm + tri + roughness + pr_urban + nest_dist + vr_dist, data = bv3, family = binomial()) vm4 = glm(bin ~ elev + vrm + tri + roughness + pr_urban + nest_dist + vr_dist, data = bvE, family = binomial()) vif(vm3) vif(vm4) vm5 = glm(bin ~ elev + vrm + roughness + pr_urban + nest_dist + vr_dist, data = bv3, family = binomial()) vm6 = glm(bin ~ elev + vrm + roughness + pr_urban + nest_dist + vr_dist, data = bvE, family = binomial()) vif(vm5) vif(vm6) vm7 = glm(bin ~ elev + vrm + slope + pr_urban + nest_dist + vr_dist, data = bv3, family = binomial()) vm8 = glm(bin ~ elev + vrm + slope + pr_urban + nest_dist + vr_dist, data = bvE, family = binomial()) vif(vm7) vif(vm8) rm(vm1, vm2, vm3, vm4, vm5, vm6, vm7, vm8) # Tidy # retaining one of those explanatory (slope, tri and roughness) variable will be appropriate ### (b - ii) collinearity between categorical explanatory variables ### library(lattice) # SLOPE # bwplot(bv3$slope ~ bv3$land_cover, cex=0.5, xlab="Habitat Class", ylab="Slope", main="Slope by Habitat Class", par.settings=list(box.rectangle=list(col="blue"), box.umbrella=list(col="blue"), plot.symbol=list(cex=1, col="blue")), scales=list(x=list(relation="same"), y=list(relation="same")), pch="|", cex.lab=3, cex.main=3) summary(aov(bv3$slope ~ bv3$land_cover)) # ELEVATION # bwplot(bv3$elev ~ bv3$land_cover, cex=0.5, xlab="Habitat Class", ylab="Elevation (m)", main="Slope by Habitat Class", par.settings=list(box.rectangle=list(col="blue"), box.umbrella=list(col="blue"), plot.symbol=list(cex=1, col="blue")), scales=list(x=list(relation="same"), y=list(relation="same")), pch="|", cex.lab=3, cex.main=3) summary(aov(bv3$elev ~ bv3$land_cover)) # VRM # bwplot(bv3$vrm ~ bv3$land_cover, cex=0.5, xlab="Habitat Class", ylab="VRM", main="Slope by Habitat Class", par.settings=list(box.rectangle=list(col="blue"), box.umbrella=list(col="blue"), plot.symbol=list(cex=1, col="blue")), scales=list(x=list(relation="same"), y=list(relation="same")), pch="|", cex.lab=3, cex.main=3) summary(aov(bv3$vrm ~ bv3$land_cover)) # TRI # bwplot(bv3$tri ~ bv3$land_cover, cex=0.5, xlab="Habitat Class", ylab="TRI", main="Slope by Habitat Class", par.settings=list(box.rectangle=list(col="blue"), box.umbrella=list(col="blue"), plot.symbol=list(cex=1, col="blue")), scales=list(x=list(relation="same"), y=list(relation="same")), pch="|", cex.lab=3, cex.main=3) summary(aov(bv3$tri ~ bv3$land_cover)) # VULTURE RESTAURANT DISTANCE # bwplot(bv3$vr_dist ~ bv3$provisioning_cat, cex=0.5, xlab="provisioning category", ylab="VR_DIST", main="vr_dist by provisioning category", par.settings=list(box.rectangle=list(col="blue"), box.umbrella=list(col="blue"), plot.symbol=list(cex=1, col="blue")), scales=list(x=list(relation="same"), y=list(relation="same")), pch="|", cex.lab=3, cex.main=3) summary(aov(bv3$vr_dist ~ bv3$land_cover)) ### (c - i) relationships between the response vs. Continuous Variables ### # SLOPE plot(jitter(bin) ~ jitter(slope), data = bv3, las = 1, xlab = "Slope", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) plot(jitter(bin) ~ jitter(slope), data = bvE, las = 1, xlab = "Slope", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) # ELEVATION plot(jitter(bin) ~ jitter(elev), data = bv3, las = 1, xlab = "Elevation", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) plot(jitter(bin) ~ jitter(elev), data = bvE, las = 1, xlab = "Elevation", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) # VRM plot(jitter(bin) ~ jitter(vrm), data = bv3, las = 1, xlab = "VRM", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) plot(jitter(bin) ~ jitter(vrm), data = bvE, las = 1, xlab = "VRM", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) # PROPORTION URBAN plot(jitter(bin) ~ jitter(pr_urban), data = bv3, las = 1, xlab = "proportion urban", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) plot(jitter(bin) ~ jitter(pr_urban), data = bvE, las = 1, xlab = "proportion urban", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) # VULTURE RESTAURANT DISTANCE plot(jitter(bin) ~ jitter(vr_dist), data = bv3, las = 1, xlab = "minimum distance from vulture restaurant", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) plot(jitter(bin) ~ jitter(vr_dist), data = bvE, las = 1, xlab = "minimum distance from vulture restaurant", pch = 20, col = "firebrick3", ylab = "Presence/Absence", cex.lab = 1.5) ### (c - ii) relationships between the response vs. Categorical Variables ### count = table(bv3$bin, bv3$land_cover) prop.table(count,2) barplot(table(bv3$bin, bv3$land_cover), ylab="Count", main="presence vs. absence Points by Habitat Type", col=c("blue", "red"), legend = rownames(count), beside=TRUE) require(plyr) ddply(bv_3, c("bin","land_cover"), summarise, n=length(land_cover)) count2 = table(bvE$bin, bvE$land_cover) prop.table(count2,2) barplot(table(bvE$bin, bvE$land_cover), ylab="Count", main="presence vs. absence Points by Habitat Type", col=c("blue", "red"), legend = rownames(count2), beside=TRUE) require(plyr) ddply(bvE, c("bin","land_cover"), summarise, n=length(land_cover)) # It looks like bare-land area is not represented in fixes this might cause problem while running GLMM # One way to tackle it to convert the bare-land into cultivated according to the following reason # Bare-land = Lands with exposed soil, sand, or rocks and never has more than 10 % vegetated cover during any time of the year # Cultivated = Lands covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems). rm(count, count2) # Tidy ### (d) Interactions Between Explanatory Variables ### plot(bv3$pr_urban ~ bv3$elev) plot(bv3$pr_urban ~ bv3$slope) plot(bv3$pr_urban ~ bv3$vrm) #### Tidying land cover and provisioning category data before running the GLMM #### # LAND COVER # # 3 time pseudo points bv3$land_cover = as.character(bv3$land_cover) # Converting to character bv3$land_cover[bv3$land_cover %in% "Shr"] = "Gra" bv3$land_cover[bv3$land_cover %in% c("Wet", "Wat", "Urb", "Bar")] = "Oth" # Equal pseudo points bvE$land_cover = as.character(bvE$land_cover) # Converting to character bvE$land_cover[bvE$land_cover %in% "Shr"] = "For" bvE$land_cover[bvE$land_cover %in% c("Wet", "Wat", "Urb", "Bar")] = "Other" bv3$land_cover = as.factor(bv3$land_cover) # Revert back to factor bvE$land_cover = as.factor(bvE$land_cover) summary(bv3$land_cover) # Should now have 5 categories of land cover summary(bvE$land_cover) # checking conversion # require(plyr) ddply(bv3, c("bin","land_cover"), summarise, n=length(land_cover)) ddply(bvE, c("bin","land_cover"), summarise, n=length(land_cover)) #### Checking class assigned to variables #### #### str(bv3) str(bvE) #### Model building explanation #### ## Group the explanatory variables into functional groups ## # (i) terrain (Slope, elevation, vector ruggedness index) # (ii) behavioral (minimum distance from nest, minimum distance from vulture restaurant, VR provisioning category) # (iii) landscape scale habitat type (Land cover) # (iv) landscape scale urbanization (Proportion urban) # There is good evidence from other studies that terrain, behavioral parameters # are a key determinant of the use of landscape by bearded vulture. Therefore, I # included these parameters in all the models #### Hypothesis 1 what land classes is selected and avoided #### # Data set 3:1 # ML1 = glmer(bin ~ scale(elev) + scale(slope) + scale(vrm) + scale(nest_dist) + scale(vr_dist) + land_cover + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") # Global model summary(ML1) result.ML1 = dredge(ML1) # this will run all kind of combination # Data set 1:1 # ML2 = glmer(bin ~ scale(elev) + scale(slope) + scale(vrm) + scale(nest_dist) + scale(vr_dist) + land_cover + (1|Name), data = bvE, family = binomial (link = "logit"), na.action = "na.fail") # Global model summary(ML2) result.ML2 = dredge(ML2) ### Visualization of output with appropriate scale ### v1 = predictorEffects(ML1, ~ land_cover) v2 = predictorEffects(ML2, ~ land_cover) library(lattice) plot(v1, lines=list(multiline=F), grid=TRUE, lattice=list(key.args=list(space="right", border=TRUE)), main = "Land cover selection plot 3:1", axes=list(y=list(transform=exp, lab="use"))) plot(v2, lines=list(multiline=F), grid=TRUE, lattice=list(key.args=list(space="right", border=TRUE)), main = "Land cover selection plot 1:1", axes=list(y=list(transform=exp, lab="use"))) #### Hypothesis 2 Is urban area avoided #### # Data set 3:1 # MU1 = glmer(bin ~ scale(elev) + scale(slope) + scale(vrm) + scale(nest_dist) + scale(vr_dist) + scale(pr_urban_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") # Global model result.MU1 = dredge(MU1) write.csv(result.MU1,"result_dredge_pr_urban_x3.csv", row.names = F, col.names = T) # Data set 1:1 # MU2 = glmer(bin ~ scale(elev) + scale(slope) + scale(vrm) + scale(nest_dist) + scale(vr_dist) + scale(pr_urban_210) + (1|Name), data = bvE, family = binomial (), na.action = "na.fail") # Global model summary(MU2) result.MU2 = dredge(MU2) ### Visualization of output with appropriate scale for urban ### t = seq(from = 0, to = 100, length = 200) urb_3 = as.data.frame(effect("scale(pr_urban_210)", MU1, xlevels=list(pr_urban_210=t))) urb_e = as.data.frame(effect("scale(pr_urban_210)", MU2, xlevels=list(pr_urban_210=t))) plot(urb_3$pr_urban_210, urb_3$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0,0.01), main = "Proportion urban", xlab="pr urban", ylab="Use") # Add confidence interval polygon(c(urb_3$pr_urban_210, rev(urb_3$pr_urban_210)), c(urb_3$lower, rev(urb_3$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) plot(urb_e$pr_urban_210, urb_e$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0,0.4), main = "Proportion urban (b)", xlab="pr urban", ylab="Use") # Add confidence interval polygon(c(urb_e$pr_urban_210, rev(urb_e$pr_urban_210)), c(urb_e$lower, rev(urb_e$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) # According to 50% (this is taken as being avoided) I get the following pr_urban # 3:1 (0.008856799; Half fit = 0.00443; pr_urban = 33.66) # 1:1 (fit = 0.16238; pr_urban = 37.18) # I will filter the urban cover layer omitting value below these threshold then # calculate the minimum distance rm(urb_3, urb_e, t) # Tidy ### calculate minimum distance of fix/pseudo absence with urban ### urb_210 = raster("urban_mean_210.tif") # uploading layer val_210 = unique(urb_210) val_210 rcl_a_3 = matrix(c(val_210, rep(NA, 17), rep(1, 33)), ncol = 2) rcl_a_3 rcl_a_e = matrix(c(val_210, rep(NA, 34), rep(1, 65)), ncol = 2) rcl_a_e urb_a3 = reclassify(urb_210, rcl = rcl_a_3) plot(urb_a3) writeRaster(urb_a3, "reclass_urban_210_3.tif") urb_a1 = reclassify(urb_210, rcl = rcl_a_e) plot(urb_a1) writeRaster(urb_a1, "reclass_urban_300_ae.tif") rm(rcl_a_3, rcl_a_e, val_210, urb_210, urb_3) # Tidy ## Projection ## prj = "+proj=longlat +datum=WGS84 +zone=35s" utm = "+proj=utm +datum=WGS84 +zone=35s" ## Read in track data bv_3 = bv3 coordinates(bv_3) =~ X + Y # Turn into spatial coordinates crs(bv_3)= crs(prj) # Assigning projection bv_3 = spTransform(bv_3, crs(utm)) bv_E = bvE coordinates(bv_E) =~ X + Y # Turn into spatial coordinates crs(bv_E)= crs(prj) # Assigning projection bv_E = spTransform(bv_E, crs(utm)) ## Turning urban pixel into spatial points up_a3 = rasterToPoints(urb_a3, spatial = FALSE) head(up_a3) up_a3 = as.data.frame(up_a3) write.csv(up_a3,"urban_point_210_3.csv", quote = F, row.names = F , col.names = T) up_a1 = rasterToPoints(urb_a1, spatial = FALSE) head(up_a1) up_a1 = as.data.frame(up_a1) write.csv(up_a1,"urban_point_300_ae.csv", quote = F, row.names = F , col.names = T) rm(urb_a1, urb_a3, up_a1, up_a3) # Tidy ## Converting into UTM to calculate distance ## up_a3 = read.csv("urban_point_210_3.csv", header = T) coordinates(up_a3) = ~ x + y crs(up_a3) = crs(prj) up_a3 = spTransform(up_a3, crs(utm)) up_a1 = read.csv("urban_point_300_ae.csv", header = T) coordinates(up_a1) = ~ x + y crs(up_a1) = crs(prj) up_a1 = spTransform(up_a1, crs(utm)) ## Calculate distance ## urb_dist_a3 = get.knnx(data = coordinates(up_a3), query = coordinates(bv_3), k = 1, algorithm = "kd_tree") bv3$urb_dist_210 = urb_dist_a3$nn.dist/1000 urb_dist_a1 = get.knnx(data = coordinates(up_a1), query = coordinates(bv_E), k = 1, algorithm = "kd_tree") bvE$urb_dist_210 = urb_dist_a1$nn.dist/1000 rm(prj, utm, up_a3, up_a1, urb_dist_a3, urb_dist_a1, bv_3, bv_E) # Tidy write.csv(bv3,"extended_result_x3.csv", row.names = F, col.names = T) ### Checking correlation ### cor3 = cor(bv3[,c(7, 8, 9, 10, 11, 12, 14, 15, 16)], method = "pearson") round(cor3, 2) ggcorrplot(cor3, hc.order = TRUE, type = "lower", lab = TRUE) corE = cor(bvE[,c(8, 9, 10, 11, 12, 14, 16, 18)], method = "pearson") round(corE, 2) ggcorrplot(corE, hc.order = TRUE, type = "lower", lab = TRUE) rm(cor3, corE) # Tidy #### Model with distance to urban cover #### MD1 = glmer(bin ~ scale(elev) + scale(slope) + scale(vrm) + scale(nest_dist) + scale(vr_dist) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") # Global model result.MD1 = dredge(MD1) result.MD1 write.csv(result.MD1,"result_dredge_distance_urban_x3.csv", row.names = F, col.names = T) vif(MD1) # Result elevation is moderately correlated with urban distance ## step selection method ## MD2 = glmer(bin ~ scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD3 = glmer(bin ~ scale(nest_dist) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD4 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD5 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(slope) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD6 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(slope) + scale(elev) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD7 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(slope) + scale(vrm) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD8 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(elev) + scale(vrm) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD9 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(vrm) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") MD10 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(elev) + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") sel = model.sel(MD2, MD3, MD4, MD5, MD6, MD7, MD8, MD9, MD10) sel # Out put #Model selection table # (Int) scl(urb_dst_210) scl(nst_dst) scl(vr_dst) scl(slp) scl(elv) scl(vrm) family df logLik AICc delta Weight #MD7 -4.752 0.07973 -4.780 -0.7489 0.6781 0.5954 binomial(logit) 7 -31452.46 62918.9 0.00 1 #MD6 -4.902 -0.05551 -4.743 -0.9719 0.7985 0.3266 binomial(logit) 7 -32261.12 64536.2 1617.33 0 #MD5 -4.925 0.09401 -4.931 -0.8191 0.8293 binomial(logit) 6 -32452.08 64916.2 1997.25 0 #MD8 -4.796 -0.04045 -4.828 -0.8614 0.3452 0.8246 binomial(logit) 7 -32537.23 65088.5 2169.54 0 summary(MD7) plot(allEffects(MD7)) rm(MD2, MD3, MD4, MD5, MD6, MD7, MD8, MD9, MD10) # Tidy #### final model with land cover classes #### M1 = glmer(bin ~ scale(nest_dist) + scale(vr_dist) + scale(slope) + scale(vrm) + land_cover + scale(urb_dist_210) + (1|Name), data = bv3, family = binomial (link = "logit"), na.action = "na.fail") summary(M1) ### Graphical parameter ### par(pty='s') # visualizing slope # ts = seq(from = 0, to = 80, length = 80) s = as.data.frame(effect("scale(slope)", M1, xlevels=list(slope = ts))) plot(s$slope, s$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0.003,0.25), xlab="Slope (degree)", ylab="Probability of use") # Add confidence interval polygon(c(s$slope, rev(s$slope)), c(s$lower, rev(s$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) box() # visualizing vector ruggedness index # tr = seq(from = 0, to = 0.4, length = 50) r = as.data.frame(effect("scale(vrm)", M1, xlevels=list(vrm = tr))) plot(r$vrm, r$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0,1), xlab="vrm", ylab="Probability of use") # Add confidence interval polygon(c(r$vrm, rev(r$vrm)), c(r$lower, rev(r$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) box() # Visualizing nest distance # tn = seq(from = 0, to = 50, length = 50) nest = as.data.frame(effect("scale(nest_dist)", M1, xlevels=list(nest_dist = tn))) plot(nest$nest_dist, nest$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0,1), xlab="Distance from nest (km)", ylab="Probability of use") # Add confidence interval polygon(c(nest$nest_dist, rev(nest$nest_dist)), c(nest$lower, rev(nest$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) box() # visualizing vulture restaurant distance # tv = seq(from = 0, to = 100, length = 100) vr = as.data.frame(effect("scale(vr_dist)", M1, xlevels=list(vr_dist = tv))) plot(vr$vr_dist, vr$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0.001,0.035), xlab="Distance from vulture restaurant (km)", ylab="Probability of use") # Add confidence interval polygon(c(vr$vr_dist, rev(vr$vr_dist)), c(vr$lower, rev(vr$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) box() # visualizing urban distance # tu = seq(from = 0, to = 34, length = 34) u = as.data.frame(effect("scale(urb_dist_210)", M1, xlevels=list(urb_dist_210 = tu))) plot(u$urb_dist_210, u$fit, col="black", cex=0.5, las=1, type="l", bty="n", ylim=c(0.004,0.02), xlab="Distance from built-up area (km)", ylab="Probability of use") # Add confidence interval polygon(c(u$urb_dist_210, rev(u$urb_dist_210)), c(u$lower, rev(u$upper)), col = adjustcolor(c("black"), alpha.f = 0.2), border = NA, pch=0.2) box() # Visualizing land cover # par(pty='s') l = predictorEffects(M1, ~ land_cover) l = as.data.frame(l) l$order = c(1:4) plot(l$order, l$land_cover.fit, las=1,bty="n", pch = 16, ylim = c(0.001, 0.014), xaxt = "n", xlab = "", ylab="Probability of use") install.packages("plotrix") require(plotrix) plotCI(l$order, l$land_cover.fit, ui= l$land_cover.upper, li=l$land_cover.lower, pch = 16, xaxt = "n", xlab = "", ylab="Probability of use") segments(l$order, l$land_cover.lower, l$order, l$land_cover.upper) axis(side=1, at = c(1:4), labels = c("cultivated", "forest", "grassland", "other"), las=2) #### computing Area Under Curve #### prob = predict(M1, type=c("response")) pred = prediction(prob, bv3$bin) perf = performance(pred, measure = "tpr", x.measure = "fpr") par(pty='s') # this makes a square plot plot(perf, colorize = TRUE, main="ROC curve Admissions", xlab="False positive", ylab="True positive") abline(0, 1) text(x = 0.2, y = 0.8, labels = "AUC = 0.9868") # calculate area under curve # library(pROC) preds=predict(M1, bv3, allow.new.levels = TRUE, type="response") roc=pROC::roc(bin ~ preds, data = bv3, quiet=T) auc(roc) library(coefplot) coefplot(M1, intercept=F, main="Fixed effect coefficient") coefplot(M2, intercept=F, main="Fixed effect coefficient")