require(plyr) require(rworldmap) require(rgeos) require(fields) require(foreach) require(parallel) require(letsR) require(rgdal) require(raster) require(parallel) require(doParallel) require(foreach) require(rgeos) require(classInt) library(gstat) library(RColorBrewer) require(maps) require(car) ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Summarize BBS Data to Species x Route ############################################################################################################################## ###### Calculate mean abundances at locations # Load in curated BBS Dataset load('BBS.Rdata') # We start with a BBS data set where the following have been added to the original data: # A unique Route ID based on the combination of countrynum, statenum, and Route # A new Species id following the included SppList file # The SpeciesTotal column from the original BBS data is used as abundance at a route # lati and longi are also from original data # Load in a table of Unique Routes 'AllUniqueRoutes' # and table of list of their surveyed years 'SurveyedYearsAll' load('BBSYearCounts.RData') # AllUnique Routes is compiled from original BBS data and includes # countrynum, statenum, Route # Our new RouteID and added YearCount # Restrict BBS to years from 1984 to 2014 BBS <- BBS[which(BBS$Year > 1984), ] # Extract unique Species IDs USp <- unique(BBS$SppID) # Initiate main data table - BBS summarized for each species x route combo Data <- as.data.frame(matrix(NA, nrow = 1, ncol = 8)) names(Data) <- c('SppID', 'RouteID', 'lat', 'long', 'Abun', 'YearCount') # Loop through species for (i in 1:length(USp)) { # Extract current species current <- BBS[which(BBS$SppID == USp[i]), ] # Which routes does it occur at? URt <- unique(current$RouteID) # Loop through routes for (j in 1:length(URt)) { # initiate chunk for species x route combo curDat <- as.data.frame(matrix(NA, nrow = 1, ncol = 8)) names(curDat) <- c('SppID', 'RouteID', 'lat', 'long', 'Abun', 'SD', 'CV', 'YearCount') # Get current route curRt <- current[which(current$RouteID == URt[j]), ] # Add Species ID curDat$SppID <- curRt$SppID[1] # Add route id curDat$RouteID <- curRt$RouteID[1] # add route lat and lon curDat$lat <- curRt$Lati[1] curDat$long <- curRt$Longi[1] # Calculate average abundance at the route # Zero abundances are added for years the route was surveyed, # but the species wasn't observed Abuns <- c(curRt$SpeciesTotal, rep(0, length(which(SurveyedYearsAll[[URt[j]]] %in% curRt$Year == F & SurveyedYearsAll[[URt[j]]] > 1984)))) curDat$Abun <- mean(Abuns) # Add number of years species was observed curDat$YearCount <- nrow(curRt) # Add to current Spp x Route to main data table Data <- rbind(Data, curDat) } print(paste0(i, ' spps dropped <(;_:)>')) } Data <- Data[2:nrow(Data), ] rm(BBS) ############################# ## Add number of years each route was surveyed Data$YearsSurveyed <- NA for (i in 1:nrow(AllUniqueRoutes)) { CurYears <- SurveyedYearsAll[[i]] CurYears <- CurYears[which(CurYears > 1984)] Data$YearsSurveyed[which(Data$RouteID == AllUniqueRoutes$RouteID[i])] <- length(CurYears) } ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Get geographic centers and calculate distances to center ############################################################################################################################## # Load in species list (file names for map shapefiles) load('SppList.Rdata') # Get unique species UniSpp <- unique(BBS$SppID) # initiate data table of geogrpahic midpoints for each species Mids <- as.data.frame(matrix(NA, nrow = length(UniSpp), ncol = 7)) names(Mid) <- c('SppID', 'x.PC', 'y.PC', 'x.GM', 'y.GM', 'x.CMD', 'y.CMD') Mids$SppID <- UniSpp # Initiate columns to store distances in the main Data table Data$DistCMD <- Data$DistGM <- Data$DistPC <- NA # Loop through species for (i in nrow(Mids)) { # get current species shapefile file name name <- SppList$scinamemap[which(SppList$SppID == Mids$SppID[i])][1] # Check if it has a range map file if (file.exists(paste0('Bird Ranges/', name, '.shp'))) { # load BirdLife international shapefile (all species stored in one folder) SppMap <- readOGR('Bird Ranges/', layer = name, verbose = F) # Restrict to known breeding range SppMap <- SppMap[SppMap$PRESENCE != 3 & SppMap$PRESENCE != 4 & SppMap$PRESENCE != 5 & SppMap$PRESENCE != 6 & SppMap$SEASONAL != 3 & SppMap$SEASONAL != 5 & SppMap$SEASONAL != 4, ] # create Presence absence matrix for use in letsR curPreAb <- lets.presab(SppMap) # calculate midpoints using three methods # PC = geographic centroid # GM = the midpoint of latitude and longtitude # CMD = the point within the range minimizing distance to all other points curMid <- lets.midpoint(curPreAb, method = 'PC') Mids$x.PC <- curMid$x Mids$y.PC <- curMid$y curMid <- lets.midpoint(curPreAb, method = 'GM') Mids$x.GM <- curMid$x Mids$y.GM <- curMid$y curMid <- lets.midpoint(curPreAb, method = 'CMD') Mids$x.CMD <- curMid$x Mids$y.CMD <- curMid$y # get BBS data for current species current <- Data[which(Data$SppID == Mids$SppID[i]), ] currentPos <- which(Data$SppID == Mids$SppID[i]) # calculate geographic distances Data$DistPC[currentPos] <- rdist.earth(current[, c('long', 'lat')], Mids[i, c('x.PC', 'y.PC')]) Data$DistGM[currentPos] <- rdist.earth(current[, c('long', 'lat')], Mids[i, c('x.GM', 'y.GM')]) Data$DistCMD[currentPos] <- rdist.earth(current[, c('long', 'lat')], Mids[i, c('x.CMD', 'y.CMD')]) } } ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Calcualte distances to the edge ############################################################################################################################## # We do not calculate distancess to edge created by lakes contained entirely in species distributions # The following 'FillTheHoles' function is required to remove these lakes # Fill the hole function # Remove holes or non-holes from a shapefile below a certain size. # hole = boolean. When TRUE will remove holes. When FALSE will remove non-holes polygons. # area = numeric. Should be one value indicating the size of the polygons above which they will be kept. FillTheHoles <- function(ranges, area = NULL, hole = TRUE) { if (is.null(area)) { area <- max(area(ranges)) + 1 } areas <- lapply(ranges@polygons, function(y){sapply(slot(y, "Polygons"), slot, "area")}) holes <- lapply(ranges@polygons, function(y){sapply(slot(y, "Polygons"), slot, "hole")}) pos <- areas N <- length(areas) keep <- logical(N) for (i in 1:N) { pos[[i]] <- (areas[[i]] < rep(area, length(areas[[i]]))) & (holes[[i]] == hole) keep[i] <- any(!pos[[i]] & !holes[[i]]) } ranges <- ranges[keep, ] N <- N - sum(!keep) pos <- pos[keep] # Remove them for (i in 1:N) { ranges@polygons[[i]]@Polygons <- ranges@polygons[[i]]@Polygons[!pos[[i]]] } while (any(sapply(ranges@polygons, function(x){x@Polygons[[1]]@hole}))) { for (i in 1:length(ranges@polygons)) { if (ranges@polygons[[i]]@Polygons[[1]]@hole) { ranges@polygons[[i]]@Polygons <- ranges@polygons[[i]]@Polygons[-1] } } } comment(ranges) <- rep("This polygon was made by BV", length(ranges)) # Remove internal borders result <- rgeos::gUnaryUnion(ranges) # Voilá! return(result) } ############################### ## Create Border Rasters and extract coordinates of edge cells ## Next we will rasterize the border of each specie's range map # This takes some time, so this process has been parallelized cl <- makeCluster(6, type = 'PSOCK', outfile = 'debug.txt', verbose = T) # feed in necessary packages clusterEvalQ(cl, require(rgdal)) clusterEvalQ(cl, require(rworldmap)) clusterEvalQ(cl, require(raster)) clusterEvalQ(cl, require(geosphere)) clusterEvalQ(cl, library(gstat)) clusterEvalQ(cl, library(rgeos)) clusterEvalQ(cl, require(stringr)) # Remember to save the FillTheHoles function above ^ clusterEvalQ(cl, source('Fill_the_hole.R')) # This is the function that will creatae the border rasters BorderMap <- function(i) { # load lakes and land map for plugging holes # These can be downloaded from https://www.naturalearthdata.com/downloads/ lakes <- readOGR('Landscape features/ne_10m_lakes', layer = 'ne_10m_lakes', verbose = F) land <- readOGR('Landscape features/ne_10m_land', layer = 'ne_10m_land', verbose = F) ##Trim to remove small lakes lakes <- FillTheHoles(lakes, area = .1, hole = F) LakesBuf <- gBuffer(lakes, width = .3) ## The St Lawrence river is included in many range maps, creating internal edges ## Make Plug to fill this gap Coords <- matrix(c(-67.72, 49.67, -66.91, 48.93, -76.15, 43.94, -76.76, 44.41, -67.72, 49.67), ncol = 2, byrow = T) bb <- Polygon(Coords) plug <- SpatialPolygons(list(Polygons(list(bb), ID = 'a')), proj4string = CRS(proj4string(land))) # Load in species list and get list of range map file names load('SppList.Rdata') mydata <- unique(SppList$scinamemap) refras <- raster(xmn = -179.8241, xmx = 180.1759, ymn = -78.44269, ymx = 79.55731, resolution = 0.04166667) ### Load species map if (file.exists(paste0('Bird Ranges/', mydata[i], '.shp'))) { CurrentSpp <- readOGR(dsn = 'Bird Ranges/', layer = paste(mydata[i]), verbose = F) # Restrict to known breeding distribution CurrentMap <-CurrentSpp[CurrentSpp$PRESENCE != 3 & CurrentSpp$PRESENCE != 4 & CurrentSpp$PRESENCE != 5 & CurrentSpp$PRESENCE != 6 & CurrentSpp$SEASONAL != 3 & CurrentSpp$SEASONAL != 5 & CurrentSpp$SEASONAL != 4, ] # Make sure there' something left if (nrow(CurrentMap) != 0) { # some maps are broken and will error out with unionSpatialPolygons test <- tryCatch(unionSpatialPolygons(CurrentMap,ID = rep(1, times = length(CurrentMap)), avoidGEOS = F), error=function(cond) { return(1) }) # make sure that didn't happen if (is.numeric(test) == T) { } else { CurrentMap <- test ## Plug the st lawrence river! plugT <- gIntersection(plug, CurrentMap, byid = T) # this code checks where and how the range map intersects with # the plug and avoids adding extra area beyond the original range map if (!is.null(plugT)) { if (length(plugT@polygons[[1]]@Polygons) == 1) { clist <- plugT@polygons[[1]]@Polygons[[1]]@coords } else { clist <- do.call(rbind, sapply(plugT@polygons[[1]]@Polygons, slot, 'coords')) } Top <- clist[which(clist[, 2] == max(clist[, 2])), ] if (class(Top)[1] == 'matrix') {Top <- Top[1, ]} Right <- clist[which(clist[, 1] == max(clist[, 1])), ] if (class(Right)[1] == 'matrix') {Right <- Right[1, ]} Bot <- clist[which(clist[, 2] == min(clist[, 2])), ] if (class(Bot)[1] == 'matrix') {Bot <- Bot[1, ]} Left <- clist[which(clist[, 1] == min(clist[, 1])), ] if (class(Left)[1] == 'matrix') {Left <- Left[1, ]} Coords <- matrix(c(Top, Right, Bot, Left, Top), ncol = 2, byrow = T) bb <- Polygon(Coords) Newplug <- SpatialPolygons(list(Polygons(list(bb), ID = 'a')), proj4string = CRS(proj4string(land))) CurrentMap <- rbind(CurrentMap, Newplug) CurrentMap <- gBuffer(CurrentMap, width = 0) CurrentMap <- gSimplify(CurrentMap, tol = 0.000001) } # Fill lakes and small holes Holess <- FillTheHoles(CurrentMap, area = NULL, hole = T) # uses lakes from earlier inlakes <- gIntersection(LakesBuf, Holess) if (!is.null(inlakes)) { JShape <- aggregate(rbind(CurrentMap, inlakes, makeUniqueIDs = T)) CurrentMap <- FillTheHoles(JShape, 0.1, hole = T) } else { CurrentMap <- FillTheHoles(CurrentMap, 0.1, hole = T) } # Separate different polygons (disjunct pieces of range) into different units # this can cause errors with some maps test <- tryCatch(disaggregate(CurrentMap), error=function(cond) { return(1) }) # make sure that didn't happen if (is.numeric(test) == T) { } else { CurrentMap <- test # Rasterize range map raz <- rasterize(CurrentMap, refras) # Creates a raster of the edge cells bounds <- boundaries(raz) # extract the latitude and longitude from edge cells Edges <- as.data.frame(xyFromCell(bounds, cell = which(values(bounds) == 1))) Edges$Poly <- values(raz)[which(values(bounds) == 1)] Edges$Cell <- which(values(bounds) == 1) # Get info on interior cells as well Inners <- as.data.frame(xyFromCell(bounds, cell = which(values(bounds) == 0))) Inners$Poly <- values(raz)[which(values(bounds) == 0)] Inners$Cell <- which(values(bounds) == 0) # save everything for later use eval(parse(text = paste0("save(raz, bounds, Edges, Inners, file = 'BorderRasters/", mydata[i], "_BorderRast.Rdata')"))) } } } } rm(list = setdiff(ls(), lsf.str())) } # Run border map function in cluster across species # there are 430 species with range map files i <- c(1:430) clusterApplyLB(cl, i, BorderMap) ############################### ## Use outputs from the BorderMap function ## to calculate from BBS location to range edges # Initiate column in main BBS data table for edge distances Data$DistEdge <- NA # Get list of unique species UniqueSpp <- unique(Data$SppID) for (i in 1:length(UniqueSpp)) { # get position observations of current species and extract them curPos <- which(Data$SppID == UniqueSpp[i]) CurrentSpp <- Data[curPos, ] # make sure border raster was successful for the current species if (file.exists(paste0('BorderRasters/', SppList$scinamemap[which(SppList$SppID == UniqueSpp[i])][1], '_BorderRast.Rdata'))) { # open border raster for current species load(paste0('BorderRasters/', SppList$scinamemap[which(SppList$SppID == UniqueSpp[i])][1], '_BorderRast.Rdata')) # Check to see if BBS locations are within or outside the birdlife rangemap Ins <- extract(raz, CurrentSpp[, c('long', 'lat')]) # Loop through routes for (j in 1:nrow(CurrentSpp)) { # if Ins is NA, the location is outside range # calculate distance to all edge cells and choose the minimum # make the distance negative for outside location if(is.na(Ins[j])) { CurrentSpp$DistEdge[j] <- -1*min(distGeo(CurrentSpp[j, 4:3], Edges[, 1:2])) } else { # For BBS locations within the range, we can reduce computations time # by only calculating distance for the edges of the polygon # where the location falls curEdges <- Edges[which(Edges$Poly == Ins[j]), ] CurrentSpp$DistEdge[j] <- min(distGeo(CurrentSpp[j, c('long', 'lat')], curEdges[, 1:2])) } } # Fill in edge distances for current species Data$DistEdge[curPos] <- CurrentSpp$DistEdge } print(paste(i, 'spec-locs dist-racted')) } ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Generate occupied climate space and calculate climate position for BBS locations ############################################################################################################################## # In climate space, we calculate euclidean distances using the following function # this code can be save so that it can be run in parallel # Generate euclidean distance function euclidean <- function(p1, p2) { if (!is.matrix(p1)) { p1 <- matrix(p1, ncol = 2) } if (!is.matrix(p2)) { p2 <- matrix(p2, ncol = 2) } lp1 <- nrow(p1) lp2 <- nrow(p2) # if different size increase the size of the smaller if (lp1 != lp2) { if (lp1 > lp2) { if (lp1 %% lp2 != 0) { stop("nrow(p2) is not a multiple of nrow(p1)") } else { p2 <- do.call(rbind, replicate((lp1 / lp2), p2, simplify=FALSE)) } } if (lp1 < lp2) { if (lp2 %% lp1 != 0) { stop("nrow(p1) is not a multiple of nrow(p2)") } else { p1 <- do.call(rbind, replicate((lp2 / lp1), p1, simplify=FALSE)) } } } lp1 <- nrow(p1) dista <- numeric(lp1) # Calculate distance for (i in 1:lp1) { dista[i] <- as.vector(dist(rbind(p1[i, ], p2[i, ]))) } return(dista) } # This function will be used to calculate distances to the edge of occupied # climate space # code should be saved so that it can be run in parralel distBord <- function(r3, current, showPlot = TRUE) { # r3 = raster with enviromental values of presence (1) or absence (0) # current = data.frame with the information on the coordinates of the climate axes # named PC1 and PC2, and with abundance values named Abun. # showPlot if True will generate the plots, if false nothing is ploted. # Create boundary raster rNA <- r3 values(rNA)[which(values(rNA == 0))] <- NA boundRast <- boundaries(rNA, type = 'inner') # Get coordinates coords <- cbind(current$PC1, current$PC2) dist_line <- NA for(i in 1:nrow(current)) { if (!is.na(extract(boundRast, matrix(coords[i, ], ncol = 2))) & extract(boundRast, matrix(coords[i, ], ncol = 2)) == 1) { dist_line[i] <- 0 } else { dist_line[i] <- min(euclidean(coords[i, ], xyFromCell(boundRast, which(values(boundRast) == 1))), na.rm = T) if (is.na(extract(boundRast, matrix(coords[i, ], ncol = 2)))) { dist_line[i] <- dist_line[i]*-1 } } } # Plot dist_bord <- dist_line ii <- cut(log10(current$Abun), breaks = seq(min(log10(current$Abun)), max(log10(current$Abun)), len = 100), include.lowest = TRUE) ## Use bin indices, ii, to select color from vector of n-1 equally spaced colors if (showPlot) { par(mfrow = c(1, 2)) colors <- colorRampPalette(c("blue", "red"))(99)[ii] plot(r3, legend = F, xlab = "PC1", ylab = "PC2") plot(boundRast, add = T) points(coords, col = adjustcolor(colors, .4), pch = 20) plot(y = current$Abun, x = dist_bord, xlab = "Distance to border", ylab = "Abundance") abline(lm(current$Abun ~ dist_bord)) } return(dist_line) } # Calculations in climate space for all species can take some time # so we run them in parallel # Set up cluster cl <- makeCluster(6, type = 'PSOCK', outfile = 'debug.txt', verbose = T) # Push resources out to cluster clusterEvalQ(cl, library(rgeos)) clusterEvalQ(cl, library(geosphere)) clusterEvalQ(cl, library(rgdal)) clusterEvalQ(cl, library(raster)) clusterEvalQ(cl, library(letsR)) # Fill the hole function will be needed - see previous section clusterEvalQ(cl, source('Fill_the_hole.R')) # also the euclidean and distbord functions clusterEvalQ(cl, source('euclidean.R')) clusterEvalQ(cl, source('distbord.R')) # this is the main function that generates occupied cliamte space # calculates centers and climate distances from BBS locations to centers and edge ####################### ClimDistEm <- function(i) { ################################## # This is a raster stack of worldclim 2.0 variables # including monthly temp and precips as well as bioclim variables # can be downloaded at https://www.worldclim.org/ load('YearClimStack.Rdata') # initiate new columns to collect climate data for each BBS route x species observation ndcol <- ncol(Data) Data[, (ndcol + 1):(ndcol + 6)] <- NA names(Data)[(ndcol + 1):(ndcol + 6)] <- c('dPC', 'dGM','dCMD', 'dEdge', 'PC1', 'PC2') Spp <- unique(Data$SppID) # get data for current species curData <- Data[which(Data$SppID == Spp[i]), ] # initiate data table to collect climate data for each species SppSum <- as.data.frame(matrix(NA, nrow = 1, ncol = 15)) names(SppSum) <- c('SppID', 'Common.name', 'Order', 'Family', 'NLoc.Cent', 'PC1min', 'PC1max', 'PC2min', 'PC2max', 'PC1.PC', 'PC2.PC', 'PC1.GM', 'PC2.GM', 'PC1.CMD', 'PC2.CMD') SppSum$SppID <- Spp[i] SppSum$Common.name <- SppList$Common.Name[which(SppList$SppID == Spp[i])][1] SppSum$Order <- SppList$Order[which(SppList$SppID == Spp[i])][1] SppSum$Family <- SppList$Family[which(SppList$SppID == Spp[i])][1] # is there data for the species if (nrow(curData) >= 1) { ###################################### ## Extract from raster stack # Get map file name and open it sciname <- SppList$scinamemap[which(SppList$SppID == Spp[i])][1] mapdata <- rgdal::readOGR("Range Maps/", layer = sciname, verbose = F) # Restrict to known breeding distribution map <-mapdata[mapdata$PRESENCE != 3 & mapdata$PRESENCE != 4 & mapdata$PRESENCE != 5 & mapdata$PRESENCE != 6 & mapdata$SEASONAL != 3 & mapdata$SEASONAL != 5 & mapdata$SEASONAL != 4, ] # combine to one polygon oneshape1 <- unionSpatialPolygons(map, IDs = map$SISID) # equal area projection oneshape1 <- spTransform(oneshape1, CRS('+proj=wag4 +lon_0=0')) # convert BBS locations into spatial points Pts <- SpatialPoints(curData[, c('long', 'lat')], CRS('+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0')) Pts <- spTransform(Pts, CRS('+proj=wag4 +lon_0=0')) # extract climate data from birdlife range map Clims <- extract(YearRP, oneshape1, cellnumbers = T)[[1]] # pull out non NA values for mean annual temp (column 38) and annual precipitation (49) Clims <- na.omit(Clims[, c(1, 38, 49)]) # for the very few zero precipitation cells, assign the non-zero minimum for log transformation Clims[which(Clims[, 3] == 0), 3] <- min(Clims[which(Clims[, 3] > 0), 3]) # log transform precip Clims[, 3] <- log(Clims[, 3]) ## Extract Temp and Precip for BBS locations ClimsPts <- extract(YearRP, Pts, cellnumbers = T) ClimsPts <- ClimsPts[, c(1, 38, 49)] ClimsPts[, 3] <- log(ClimsPts[, 3]) # Center and scale these values to the range of values extracted from the range map ClimsPts[, 2] <- (ClimsPts[, 2] - mean(Clims[, 2]))/sd(Clims[, 2]) ClimsPts[, 3] <- (ClimsPts[, 3] - mean(Clims[, 3]))/sd(Clims[, 3]) # Now center and scale the values from the range map Clims[, 2] <- scale(Clims[, 2]) Clims[, 3] <- scale((Clims[, 3])) # PC1 = position climate in temperature # PC2 = position climate in precipitation curData[, c('PC1', 'PC2')] <- ClimsPts[, 2:3] ## Make reference raster of 0.2 resolution r <- raster(xmn = min(na.omit(c(Clims[, 2], ClimsPts[, 2])) - 1), xmx = max(na.omit(c(Clims[, 2], ClimsPts[, 2])) + 1), ymn = min(na.omit(c(Clims[, 3], ClimsPts[, 3])) - 1), ymx = max(na.omit(c(Clims[, 3], ClimsPts[, 3])) + 1), res = .2) values(r) <- 0 # Find out which cells overlap with the values extracted from the range cells <- extract(r, Clims[, 2:3], cellnumbers = T) # and assign those as occupied values(r)[cells[, 1]] <- 1 # Get 'coordinates' from entire raster xy <- xyFromCell(r, 1:ncell(r)) r1 <- r2 <- r values(r1) <- xy[, 1] values(r2) <- xy[, 2] ## and use them to create presence absence matrix for letsR functions pamENV <- cbind(values(r1), values(r2), values(r)) colnames(pamENV) <- c("x", "y", "sp") # Create the Presence absence object PAMenv <- list(pamENV, r, colnames(pamENV)[3]) names(PAMenv) <- c("Presence_and_Absence_Matrix", "Richness_Raster", "Species_name") class(PAMenv) <- "PresenceAbsence" ############################################################################################ ## Distance to Centers ## Get centers SppSum[1, 10:11] <- lets.midpoint(PAMenv, method = "PC", planar = T)[, -1] SppSum[1, 17:18] <- lets.midpoint(PAMenv, method = "GM", planar = T)[, -1] SppSum[1, 24:25] <- lets.midpoint(PAMenv, method = "CMD", planar = T)[, -1] ## Calculate distances # Generate euclidean distance function euclidean <- function(p1, p2) { if (!is.matrix(p1)) { p1 <- matrix(p1, ncol = 2) } if (!is.matrix(p2)) { p2 <- matrix(p2, ncol = 2) } lp1 <- nrow(p1) lp2 <- nrow(p2) # if different size increase the size of the smaller if (lp1 != lp2) { if (lp1 > lp2) { if (lp1 %% lp2 != 0) { stop("nrow(p2) is not a multiple of nrow(p1)") } else { p2 <- do.call(rbind, replicate((lp1 / lp2), p2, simplify=FALSE)) } } if (lp1 < lp2) { if (lp2 %% lp1 != 0) { stop("nrow(p1) is not a multiple of nrow(p2)") } else { p1 <- do.call(rbind, replicate((lp2 / lp1), p1, simplify=FALSE)) } } } lp1 <- nrow(p1) dista <- numeric(lp1) # Calculate distance for (i in 1:lp1) { dista[i] <- as.vector(dist(rbind(p1[i, ], p2[i, ]))) } return(dista) } # Calculate distance of BBS locations from centers curData$dPC <- euclidean(as.matrix(SppSum[1, 10:11]), as.matrix(curData[, 16:17])) curData$dGM <- euclidean(as.matrix(SppSum[1, 17:18]), as.matrix(curData[, 16:17])) curData$dCMD <- euclidean(as.matrix(SppSum[1, 24:25]), as.matrix(curData[, 16:17])) # Get the minimum and maximum temperature and precipitation values from the species distribution SppSum$PC1min <- min(Clims[, 2], na.rm = T) SppSum$PC1max <- max(Clims[, 2], na.rm = T) SppSum$PC2min <- min(Clims[, 3], na.rm = T) SppSum$PC2max <- max(Clims[, 3], na.rm = T) ## Center sample size SppSum$NLoc.Cent <- length(!is.na(curData$dCMD)) ############################################################################################ ## Distance to Edge edgeData <- curData pos <- which(!is.na(edgeData$PC1) & !is.na(edgeData$PC2)) # use the distBord function to calulate distances to edge curD <- distBord(r, edgeData[pos, ], showPlot = F) edgeData$dEdge[pos] <- curD ## sampl size edge SppSum$NLoc.Edge <- length(which(!is.na(edgeData$dEdge))) # save raster of species occupied climate range save(r, file = paste0('Climate Ranges NewRevis/ClimRangeSpp', Spp[i], '.Rdata')) # Data the data chunk for the current species save(curData, file = paste0('Climate Data By spp New Revis/ClimDataSpp', Spp[i], '.Rdata')) # save the summary data for the current species save(SppSum, file = paste0('Climate Summaries by spp NewRevis/ClimSumSpp', Spp[i], '.Rdata')) # Save the extracted climate range save(Clims, file = paste0('Spp Full Clims/FullClim_', Spp[i], '.Rdata')) } } # Run the function in parralel # There are 419 species in the data i <- 1:419 clusterApplyLB(cl, i, ClimDistEm) stopCluster(cl) ## Now data from species run in parralel needs to be stitched together ### Stich together Clim data FileList <- list.files('Climate Data By spp New Revis') load(paste0('Climate Data By spp New Revis/', FileList[1])) ClimData <- curData for (i in 2:length(FileList)) { load(paste0('Climate Data By spp New Revis/', FileList[i])) ClimData <- rbind(ClimData, curData) } ### Stich together Species summaries FileList <- list.files('Climate Summaries by spp NewRevis') load(paste0('Climate Summaries by spp NewRevis/', FileList[1])) SppSums <- SppSum for (i in 2:length(FileList)) { load(paste0('Climate Summaries by spp NewRevis/', FileList[i])) SppSums <- rbind(SppSums, SppSum) } ### Stich together Clim data FileList <- list.files('Spp Full Clims') # For converting abundance core location in scaled climate to raw climate values ClimStats <- as.data.frame(matrix(NA, nrow = length(FileList), ncol = 5)) names(ClimStats) <- c('SppID', 'MeanPC1', 'SDPC1', 'MeanPC2', 'SDPC2') for (i in 1:nrow(ClimStats)) { load(paste0('Spp Full Clims/', FileList[i])) Spliz <- strsplit(FileList[i], c('_'), fixed = T) ClimStats$SppID[i] <- as.numeric(strsplit(Spliz[[1]][2], c('.'), fixed = T)[[1]][1]) ClimStats$MeanPC1[i] <- mean(Clims[, 2], na.rm = T) ClimStats$SDPC1[i] <- sd(Clims[, 2], na.rm = T) ClimStats$MeanPC2[i] <- mean(Clims[, 3], na.rm = T) ClimStats$SDPC2[i] <- sd(Clims[, 3], na.rm = T) } ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Quantify geographic range abundance-range position relationships across species ## GRIDDED ## this includes calculation of the abundance core ############################################################################################################################## # Create reference raster # Used for grid aggregating abundance and distance values rbase <- raster(xmn = -171.7911, xmx = -45, ymn = 24, ymx = 80, resolution = 1) rbase <- projectRaster(rbase, crs = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'), res = c(100000,100000)) ## Reduce data based on minimum number of observations and years surveyed ## **** here a different minimum can be chosen for 'YearCount' (1 or 5 are the alternative cutoffs in the paper) GData <- Data[which(Data$YearCount >= 2 & Data$YearsSurveyed >= 5), ] # Initiate data frame to store outputes from geographic abundance vs distance to center or core relationships GSlopesG <- as.data.frame(matrix(NA, nrow = length(unique(GData$SppID)), ncol = 11)) names(GSlopesG) <- c('SppId', 'Common.name', 'nLocs', 'RhoPPC', 'RhoPC', 'RhoPGM', 'RhoGM', 'RhoPCMD', 'RhoCMD', 'RhoPAbun', 'RhoAbun') # FIll in species IDs GSlopesG$SppId <- unique(GData$SppID) # Remove yellow-chevroned parakeet - not needed for min 5 cutoff # This introduced species has data problems GSlopesG <- GSlopesG[-which(GSlopesG$SppId == 137), ] # Initiate data frame to store lon and late for centers of abundance AbunCentersG <- as.data.frame(matrix(NA, nrow = nrow(GSlopesG), ncol = 2)) names(AbunCentersG) <- c('lon', 'lat') # These species have problems with maps NAids <- c(132, 118, 61, 4, 3, 20, 182, 101, 153, 353, 411) for (i in 1:nrow(GSlopesG)) { # Make sure it's not a problem species if (GSlopesG$SppId[i] %in% NAids == F) { # Extract data for current species current <- GData[which(GData$SppID == GSlopesG$SppId[i]), ] # Make sure ther's data to analyze if (nrow(current) > 1 & !is.na(current$DistCMD[1])) { # convert BBS locations to spatial points SpP <- SpatialPoints(current[, c('long', 'lat')], proj4string = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) SpP <- spTransform(SpP, CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) # Use reference raster to aggregate abundances at cell level PRast <- rasterize(SpP, rbase, (current$Abun), fun = mean) # Use reference raster to aggregate distances to centers at cell level DistRast1 <- rasterize(SpP, rbase, current$DistPC, fun = mean) DistRast2 <- rasterize(SpP, rbase, current$DistGM, fun = mean) DistRast3 <- rasterize(SpP, rbase, current$DistCMD, fun = mean) # Get lat and lon of cells uniCoords <- xyFromCell(PRast, which(!is.na(values(PRast)) & !is.na(values(DistRast1)))) uniCoordsLL <- SpatialPoints(uniCoords[, 1:2], proj4string = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) uniCoordsLL <- spTransform(uniCoordsLL, CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) # extract mean aggregated abundances and distances MeanBuns <- values(PRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast1)))] MeanPC <- values(DistRast1)[which(!is.na(values(PRast)) & !is.na(values(DistRast1)))] MeanGM <- values(DistRast2)[which(!is.na(values(PRast)) & !is.na(values(DistRast2)))] MeanCMD <- values(DistRast3)[which(!is.na(values(PRast)) & !is.na(values(DistRast3)))] # Find the top abundances for estimating the abundance core # **** the value for the qunatile can be changed # **** a value of 0.9 here, values of 0.5 and 0.75 are also reported in the paper Top10 <- MeanBuns[which(MeanBuns >= quantile(MeanBuns, .9))] TopCoords <- uniCoords[which(MeanBuns >= quantile(MeanBuns, .9)), ] # Find the centroid of most abundant points if (length(Top10) >= 1) { if (length(Top10) >= 2) { centroid <- gCentroid(SpatialPoints(TopCoords[, 1:2], proj4string = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'))) } else { centroid <- gCentroid(SpatialPoints(matrix(TopCoords[1:2], nrow = 1, ncol = 2), proj4string = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'))) } # Save the location centerLL <- spTransform(centroid, CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) AbunCentersG[i, ] <- centerLL@coords } GSlopesG$Common.name[i] <- as.character(SppList$Common.Name[which(SppList$SppID == GSlopesG$SppId[i])][1]) # Calculate the PC center relationship if (nrow(uniCoords) >= 3 & !is.na(current$DistPC[1])) { GSlopesG$RhoPPC[i] <- cor.test(log10(MeanBuns), MeanPC, method = 'spearman')$p.value GSlopesG$RhoPC[i] <- cor.test(log10(MeanBuns), MeanPC, method = 'spearman')$estimate } # Calculate the GM center relationship if (nrow(uniCoords) >= 3 & !is.na(current$DistGM[1])) { GSlopesG$RhoPGM[i] <- cor.test(log10(MeanBuns), MeanGM, method = 'spearman')$p.value GSlopesG$RhoGM[i] <- cor.test(log10(MeanBuns), MeanGM, method = 'spearman')$estimate } # Calculate the CMD center relationship if (nrow(uniCoords) >= 3 & !is.na(current$DistCMD[1])) { GSlopesG$RhoPCMD[i] <- cor.test(log10(MeanBuns), MeanCMD, method = 'spearman')$p.value GSlopesG$RhoCMD[i] <- cor.test(log10(MeanBuns), MeanCMD, method = 'spearman')$estimate } # # Calculate the Abundance Core relationship if (nrow(uniCoords) >= 3 & !is.na(centerLL[1])) { # First calculate distances from BBS locations to the abundance core dists <- rdist.earth(centerLL@coords, current[, 4:3]) # Cell aggregate them as before and get teh values DistRast4 <- rasterize(SpP, rbase, dists[1, ], fun = mean) MeanCore <- values(DistRast4)[which(!is.na(values(PRast)) & !is.na(values(DistRast4)))] GSlopesG$RhoPAbun[i] <- cor.test(log10(MeanBuns), MeanCore, method = 'spearman')$p.value GSlopesG$RhoAbun[i] <- cor.test(log10(MeanBuns), MeanCore, method = 'spearman')$estimate } } } print(paste0(i, ' species plot-mapped. ')) } ######################################################################## ## Add Core coordinates to GSlopesG GSlopesG$CoreLat <- GSlopesG$CoreLon <- NA GSlopesG[, c('CoreLon', 'CoreLat')] <- AbunCentersG[, 1:2] ######################################################################## ## Additional info GSlopesG$Order <- NA GSlopesG$Family <- NA SppList$Order <- as.character(SppList$Order) SppList$Family <- as.character(SppList$Family) for (i in 1:nrow(GSlopesG)) { current <- SppList[which(SppList$SppID == GSlopesG$SppId[i]), ] if (nrow(current) > 0) { GSlopesG$Order[i] <- current$Order[1] GSlopesG$Family[i] <- current$Family[1] } } ####################################################### ## Now for the geo edge relationship ## Restrict data depending on cutoffs ## **** Change year min observations at route here with minimum YearCount EData <- Data[which(Data$YearCount >= 2 & Data$YearsSurveyed >= 5), ] ## Initiate dataframe to store edge relationships across species D2ESlopesG <- as.data.frame(matrix(NA, nrow = length(unique(EData$SppID)), ncol = 5)) names(D2ESlopesG) <- c('SppID', 'Common.name','nLocsCont', 'RhoCont', 'RhoPCont') D2ESlopesG$SppID <- unique(EData$SppID) for (i in 1:nrow(D2ESlopesG)) { # Extract Data for species current <- EData[which(EData$SppID == D2ESlopesG$SppID[i] & !is.na(EData$DistEdge)), ] D2ESlopesG$Common.name[i] <- as.character(SppList$Common.Name[which(SppList$SppID == D2ESlopesG$SppID[i])][1]) # Make sure there is data if (nrow(current) > 0) { # Convert locations into spatial points SpP <- SpatialPoints(current[, c('long', 'lat')], proj4string = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) SpP <- spTransform(SpP, CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) # Reference raster for cell aggregating PRast <- rasterize(SpP, rbase, (current$Abun), fun = mean) DistRast <- rasterize(SpP, rbase, current$DistEdge, fun = mean) # Get lat lon from cells uniCoords <- xyFromCell(PRast, which(!is.na(values(PRast)) & !is.na(values(DistRast)))) # Get the aggregated values MeanBuns <- values(PRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast)))] MeanDists <- values(DistRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast)))] D2ESlopesG$nLocsCont[i] <- nrow(uniCoords) if (nrow(uniCoords) >= 3) { D2ESlopesG$RhoCont[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate D2ESlopesG$RhoPCont[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } } # Remove yellow-chevroned parakeet, not needed for min 5 cutoff - this species causes issues D2ESlopesG <- D2ESlopesG[-which(D2ESlopesG$Common.name == 'Yellow-chevroned Parakeet'), ] ######################################################################## ## Additional info D2ESlopesG$Order <- NA D2ESlopesG$Family <- NA SppList$Order <- as.character(SppList$Order) SppList$Family <- as.character(SppList$Family) for (i in 1:nrow(D2ESlopesG)) { current <- SppList[which(SppList$SppID == D2ESlopesG$SppID[i]), ] if (nrow(current) > 0) { D2ESlopesG$Order[i] <- current$Order[1] D2ESlopesG$Family[i] <- current$Family[1] } } ######################### ## Add Core coordinates D2ESlopesG$CoreLat <- D2ESlopesG$CoreLon <- NA D2ESlopesG[, c('CoreLon', 'CoreLat')] <- AbunCentersG[, 1:2] ######################### ## Add region data ### Ecoregions ### The Nature Conservancy ecoregion maps - can be downloaded at https://geospatial.tnc.org/datasets/b1636d640ede4d6ca8f5e369f2dc368b/about EcoReg <- readOGR('terr-ecoregions-TNC/', layer = 'tnc_terr_ecoregions') EcoReg <- crop(EcoReg, c(-171.7911, -45, 24, 80)) # Consolidate into broad classifications EcoReg$Broad <- NA EcoReg$Broad[which(EcoReg$WWF_MHTNAM %in% c('Tropical and Subtropical Coniferous Forests', 'Deserts and Xeric Shrublands', 'Mediterranean Forests, Woodlands and Scrub', 'Temperate Conifer Forests'))] <- 'West' EcoReg$Broad[which(EcoReg$WWF_MHTNAM %in% c('Temperate Broadleaf and Mixed Forests', 'Tropical and Subtropical Dry Broadleaf Forests', 'Tropical and Subtropical Grasslands, Savannas and Shrublands', 'Flooded Grasslands and Savannas', 'Tropical and Subtropical Moist Broadleaf Forests') | EcoReg$ECO_NAME %in% c('Bermuda Subtropical Conifer Forests', 'Bahamian Pine Mosaic', 'East Gulf Coastal Plain', 'Florida Peninsula', 'Mid-Atlantic Coastal Plain', 'South Atlantic Coastal Plain', 'Upper West Gulf Coastal Plain', 'West Gulf Coastal Plain'))] <- 'East' EcoReg$Broad[which(EcoReg$ECO_NAME %in% c('Sierra De La Laguna Dry Forests', 'Sinaloan Dry Forests', 'Sonoran-Sinaloan Transition Subtropical Dry Forest', 'Willamette Valley - Puget Trough - Georgia Basin: Temperate Broadleaf And Mixed Forests'))] <- 'West' EcoReg$Broad[which(EcoReg$WWF_MHTNAM %in% c('Temperate Grasslands, Savannas and Shrublands') | EcoReg$ECO_NAME %in% c('Black Hills'))] <- 'Plains' EcoReg$Broad[which(EcoReg$WWF_MHTNAM %in% c('Tundra', 'Boreal Forests/Taiga'))] <- 'North' EcoReg <- spTransform(EcoReg, CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) ## Extract ecoregion data for the abundance cores CorePos <- which(!is.na(D2ESlopesG$CoreLat)) stractSP <- SpatialPoints(na.omit(D2ESlopesG[, c('CoreLon', 'CoreLat')]), CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) stractSP <- spTransform(stractSP, CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) Reg <- over(stractSP, EcoReg) # Add to species data D2ESlopesG$Reg[CorePos] <- Reg$Broad # Correct 'offshore' cores D2ESlopesG$Reg[which(1:nrow(D2ESlopesG) %in% CorePos & is.na(D2ESlopesG$Reg) & D2ESlopesG$CoreLon < -105)] <- 'West' D2ESlopesG$Reg[which(1:nrow(D2ESlopesG) %in% CorePos & is.na(D2ESlopesG$Reg) & D2ESlopesG$CoreLon > -105 & D2ESlopesG$CoreLat < 50)] <- 'East' D2ESlopesG$Reg[which(1:nrow(D2ESlopesG) %in% CorePos & is.na(D2ESlopesG$Reg) & D2ESlopesG$CoreLon > -105 & D2ESlopesG$CoreLat > 50)] <- 'North' ######################### ## Calculate distance of the core of abundance from the coast # Load in terrestrial landscape features, see above for source lakes <- readOGR('Landscape features/ne_10m_lakes', layer = 'ne_10m_lakes', verbose = F) land <- readOGR('Landscape features/ne_10m_land', layer = 'ne_10m_land', verbose = F) ##Trim to remove small lakes lakes <- FillTheHoles(lakes, area = .1, hole = F) LakesBuf <- gBuffer(lakes, width = .3) ## Make Plug - for the st lawrence river Coords <- matrix(c(-67.72, 49.67, -66.91, 48.93, -76.15, 43.94, -76.76, 44.41, -67.72, 49.67), ncol = 2, byrow = T) bb <- Polygon(Coords) plug <- SpatialPolygons(list(Polygons(list(bb), ID = 'a')), proj4string = CRS(proj4string(land))) # Fill inland water bodies WorldEdge <- gUnion(land, plug) WorldEdge <- crop(WorldEdge, extent(-180, -35, 13, 90)) WorldEdge <- FillTheHoles(WorldEdge, area = .1, hole = T) WorldEdge <- spTransform(WorldEdge, CRS("+proj=utm +zone=10 +datum=WGS84")) # Simplify the edge for faster calculations WorldBuf <- gBuffer(WorldEdge, width = -40000) WorldBuf <- spTransform(WorldBuf, CRS("+proj=longlat +datum=WGS84 +towgs84=0,0,0")) rm(WorldEdge, lakes) # initiate new columns AbunCentersG$inCont <- AbunCentersG$Dist2Cont <- NA # Calculate distances holdCore <- na.omit(AbunCentersG[, 1:2]) holdDist <- dist2Line(holdCore, WorldBuf) # Set cores outside of continent to negative distances holdover <- over(SpatialPoints(holdCore, CRS("+proj=longlat +datum=WGS84 +towgs84=0,0,0")), WorldBuf) holdover[which(is.na(holdover))] <- 0 AbunCentersG$Dist2Cont[which(!is.na(AbunCentersG$lon))] <- holdDist[, 1] AbunCentersG$inCont[which(!is.na(AbunCentersG$lon))] <- holdover AbunCentersG$Dist2Cont[which(AbunCentersG$inCont == 0)] <- AbunCentersG$Dist2Cont[which(AbunCentersG$inCont == 0)]*(-1) # Add to GSlopesG GSlopesG$CorD2Cont <- AbunCentersG$Dist2Cont GSlopesG$CorIn <- AbunCentersG$inCont # Add to D2ESlopesG D2ESlopesG$CorD2Cont <- AbunCentersG$Dist2Cont D2ESlopesG$CorIn <- AbunCentersG$inCont ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Quantify climate range abundance-range position relationships across species ## GRIDDED ## this includes calculation of the abundance core ############################################################################################################################## ## Restrict data depending on minimums ## ***** change YearCount for min counts GClimData <- ClimData[which(ClimData$YearCount >= 2 & ClimData$YearsSurveyed >= 5), ] # Initiate data frame for climate relationships SppClimsG<- as.data.frame(matrix(NA, nrow = length(unique(GClimData$SppID)), ncol = 45)) names(SppClimsG) <- c('SppID', 'Common.name', 'nLocsE', 'RhoE', 'RhoPE', 'nLocsCent', 'RhoPC', 'RhoPPC', 'SlopeGM', 'RhoGM', 'RhoPGM', 'SlopeCMD', 'RhoCMD', 'RhoPCMD', 'nLocsCore', 'RhoCore', 'RhoPCore', 'CorePC1', 'CorePC2') # Fill in species IDs SppClimsG$SppID <- unique(GClimData$SppID) # Loop through species for (i in 1:nrow(SppClimsG)) { # Get BBS data for current species current <- GClimData[which(GClimData$SppID == SppClimsG$SppID[i] & !is.na(GClimData$dEdge)), ] SppClimsG$Common.name[i] <- as.character(SppList$Common.Name[which(SppList$SppID == SppClimsG$SppID[i])][1]) # grid points by rounding lat lon coordinates current$PC1 <- round_any(current$PC1, .2) current$PC2 <- round_any(current$PC2, .2) # Get unique locations and calculate the average abundance and distance values uniCoords <- unique(current[, c('PC1', 'PC2')]) MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dEdge[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } # Number of locations SppClimsG$nLocsE[i] <- nrow(uniCoords) # Calculate edge relationships if (nrow(uniCoords) >= 3) { SppClimsG$RhoE[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClimsG$RhoPE[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } ## Centers current <- GClimData[which(GClimData$SppID == SppClimsG$SppID[i] & !is.na(GClimData$dCMD)), ] # grid points by rounding lat lon coordinates current$PC1 <- round_any(current$PC1, .2) current$PC2 <- round_any(current$PC2, .2) # get unique locations and calculate averages uniCoords <- unique(current[, c('PC1', 'PC2')]) # PC center MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dPC[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } SppClimsG$nLocsCent[i] <- nrow(uniCoords) if (nrow(uniCoords) >= 3) { SppClimsG$RhoPC[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClimsG$RhoPPC[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } # GM center MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dGM[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } if (nrow(uniCoords) >= 3) { SppClimsG$RhoGM[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClimsG$RhoPGM[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } # CMD center MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dCMD[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } if (nrow(uniCoords) >= 3) { SppClimsG$RhoCMD[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClimsG$RhoPCMD[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } ## Core current <- GClimData[which(GClimData$SppID == SppClimsG$SppID[i]), ] # grid points by rounding lat lon coordinates current$PC1 <- round_any(current$PC1, .2) current$PC2 <- round_any(current$PC2, .2) uniCoords <- unique(current[, c('PC1', 'PC2')]) # Core MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } # Get most abundance locations for calculating abundance core locations # **** change the qunatile value to test other thresholds (0.5 or 0.75) Top10 <- MeanBuns[which(MeanBuns >= quantile(MeanBuns, .9, na.rm = T))] TopCoords <- uniCoords[which(MeanBuns >= quantile(MeanBuns, .9, na.rm = T)), ] # Calculate core location and calculate distances if (length(Top10) >= 1) { centroid <- c(mean(TopCoords$PC1), mean(TopCoords$PC2)) SppClimsG[i, c('CorePC1', 'CorePC2')] <- centroid MeanDists <- sqrt((uniCoords$PC1 - centroid[1]) ^ 2 + (uniCoords$PC2 - centroid[2]) ^ 2) } else {center <- NA} # Calculate location SppClimsG$nLocsCore[i] <- nrow(uniCoords) if (any(!is.na(MeanDists))) { if (nrow(uniCoords) >= 3) { SppClimsG$RhoCore[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClimsG$RhoPCore[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } } # Add additional information SppClimsG$PC1.CMD <- SppSums$PC1.CMD SppClimsG$PC2.CMD <- SppSums$PC2.CMD SppClimsG$PC1.GM <- SppSums$PC1.GM SppClimsG$PC2.GM <- SppSums$PC2.GM SppClimsG$PC1.PC <- SppSums$PC1.PC SppClimsG$PC2.PC <- SppSums$PC2.PC SppClimsG$PC1max <- SppSums$PC1max SppClimsG$PC1min <- SppSums$PC1min SppClimsG$PC2max <- SppSums$PC2max SppClimsG$PC2min <- SppSums$PC2min ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Quantify geographic range abundance-range position relationships across species ## NON_GRIDDED ## this includes calculation of the abundance core ############################################################################################################################## ######### # See comments for the GRIDDED section for more thorough comments EData <- Data[which(Data$YearCount >= 2 & Data$YearsSurveyed >= 5), ] D2ESlopes <- as.data.frame(matrix(NA, nrow = length(unique(EData$SppID)), ncol = 4)) names(D2ESlopes) <- c('SppID', 'Common.name', 'RhoCont', 'RhoPCont') D2ESlopes$SppID <- unique(EData$SppID) for (i in 1:nrow(D2ESlopes)) { current <- EData[which(EData$SppID == D2ESlopes$SppID[i] & !is.na(EData$DistEdge)), ] D2ESlopes$Common.name[i] <- as.character(SppList$Common.Name[which(SppList$SppID == D2ESlopes$SppID[i])][1]) D2ESlopes$nLocsCont[i] <- nrow(current) if (nrow(current) >= 10) { D2ESlopes$RhoCont[i] <- cor.test(log10(current$Abun), current$DistEdge, method = 'spearman', use = 'complete')$estimate D2ESlopes$RhoPCont[i] <- cor.test(log10(current$Abun), current$DistEdge, method = 'spearman', use = 'complete')$p.value } } D2ESlopes$Order <- NA D2ESlopes$Family <- NA SppList$Order <- as.character(SppList$Order) SppList$Family <- as.character(SppList$Family) for (i in 1:nrow(D2ESlopes)) { current <- SppList[which(SppList$SppID == D2ESlopes$SppID[i]), ] if (nrow(current) > 0) { D2ESlopes$Order[i] <- current$Order[1] D2ESlopes$Family[i] <- current$Family[1] } } GData <- Data[which(Data$YearCount >= 2 & Data$YearsSurveyed >= 5), ] GSlopes <- as.data.frame(matrix(NA, nrow = length(unique(GData$SppID)), ncol = 11)) names(GSlopes) <- c('SppID', 'Common.name', 'nLocs', 'CorPC', 'WithinPC', 'CorGM', 'WithinGM', 'CorCMD', 'WithinCMD', 'CorAbun', 'WithinAbun') GSlopes$SppID <- unique(GData$SppID) AbunCenters <- as.data.frame(matrix(NA, nrow = nrow(GSlopes), ncol = 2)) names(AbunCenters) <- c('lon', 'lat') for (i in 1:nrow(GSlopes)) { if (GSlopes$SppID[i] %in% NAids == F) { current <- GData[which(GData$SppID == GSlopes$SppID[i]), ] GSlopes$nLocs[i] <- nrow(current) Top10 <- current[which(current$Abun >= quantile(current$Abun, .9)), ] if (length(Top10) >= 1) { centroid <- gCentroid(SpatialPoints(Top10[, c('long', 'lat')])) center <- centroid@coords AbunCenters[i, ] <- center } else {center <- NA} if (nrow(current) >= 10 & !is.na(current$DistPC[1])) { GSlopes$pPC[i] <- cor.test(log10(current$Abun), current$DistPC, method = 'spearman')$p.value GSlopes$CorPC[i] <- cor.test(log10(current$Abun), current$DistPC, method = 'spearman')$estimate } if (nrow(current) >= 10 & !is.na(current$DistGM[1])) { GSlopes$pGM[i] <- cor.test(log10(current$Abun), current$DistGM, method = 'spearman')$p.value GSlopes$CorGM[i] <- cor.test(log10(current$Abun), current$DistGM, method = 'spearman')$estimate } if (nrow(current) >= 10 & !is.na(current$DistCMD[1])) { GSlopes$pCMD[i] <- cor.test(log10(current$Abun), current$DistCMD, method = 'spearman')$p.value GSlopes$CorCMD[i] <- cor.test(log10(current$Abun), current$DistCMD, method = 'spearman')$estimate } if (nrow(current) >= 10 & !is.na(center[1])) { dists <- rdist.earth(center, current[, c('long', 'lat')]) GSlopes$pAbun[i] <- cor.test(log10(current$Abun), dists[1, ], method = 'spearman')$p.value GSlopes$CorAbun[i] <- cor.test(log10(current$Abun), dists[1, ], method = 'spearman')$estimate } } } ## Add Core coordinates GSlopes[, c('CoreLon', 'CoreLat')] <- AbunCenters[, 1:2] GSlopes$Order <- NA GSlopes$Family <- NA SppList$Order <- as.character(SppList$Order) SppList$Family <- as.character(SppList$Family) for (i in 1:nrow(GSlopes)) { current <- SppList[which(SppList$SppID == GSlopes$SppID[i]), ] if (nrow(current) > 0) { GSlopes$Order[i] <- current$Order[1] GSlopes$Family[i] <- current$Family[1] } } ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Quantify climate range abundance-range position relationships across species ## NON_GRIDDED ## this includes calculation of the abundance core ############################################################################################################################## ## See gridded section for more thorough comments ## vv change for min counts GClimData <- ClimData[which(ClimData$YearCount >= 2 & ClimData$YearsSurveyed >= 5), ] SppClims<- as.data.frame(matrix(NA, nrow = length(unique(GClimData$SppID)), ncol = 16)) names(SppClims) <- c('SppID', 'Common.name', 'nLocsE', 'RhoE', 'RhoPE', 'nLocsCent', 'RhoPC', 'RhoPPC', 'RhoGM', 'RhoPGM', 'RhoCMD', 'RhoPCMD', 'RhoCore', 'RhoPCore', 'CorePC1', 'CorePC2') SppClims$SppID <- unique(GClimData$SppID) for (i in 1:nrow(SppClims)) { current <- GClimData[which(GClimData$SppID == SppClims$SppID[i] & !is.na(GClimData$dEdge)), ] SppClims$Common.name[i] <- as.character(SppList$Common.Name[which(SppList$SppID == SppClims$SppID[i])][1]) uniCoords <- unique(current[, c('PC1', 'PC2')]) MeanBuns <- current$Abun MeanDists <- current$dEdge SppClims$nLocsE[i] <- length(which(!is.na(MeanDists))) if (nrow(uniCoords) >= 3) { SppClims$RhoE[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClims$RhoPE[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } ## Centers current <- GClimData[which(GClimData$SppID == SppClims$SppID[i] & !is.na(GClimData$dCMD)), ] uniCoords <- unique(current[, c('PC1', 'PC2')]) # PC MeanBuns <- current$Abun MeanDists <- current$dPC SppClims$nLocsCent[i] <- length(which(!is.na(MeanDists))) if (nrow(uniCoords) >= 3) { SppClims$RhoPC[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClims$RhoPPC[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } # GM MeanBuns <- current$Abun MeanDists <- current$dGM if (nrow(uniCoords) >= 3) { SppClims$RhoGM[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClims$RhoPGM[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } # CMD MeanBuns <- current$Abun MeanDists <- current$dCMD if (nrow(uniCoords) >= 3) { SppClims$RhoCMD[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClims$RhoPCMD[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } ## Core current <- GClimData[which(GClimData$SppID == SppClims$SppID[i]), ] uniCoords <- current[, c('PC1', 'PC2')] # Core MeanBuns <- current$Abun MeanDists <- NA Top10 <- MeanBuns[which(MeanBuns >= quantile(MeanBuns, .9, na.rm = T))] TopCoords <- current[which(MeanBuns >= quantile(MeanBuns, .9, na.rm = T)), c('PC1', 'PC2')] if (length(Top10) >= 1) { centroid <- c(mean(TopCoords$PC1), mean(TopCoords$PC2)) SppClims[i, c('CorePC1', 'CorePC2')] <- centroid MeanDists <- sqrt((uniCoords$PC1 - centroid[1]) ^ 2 + (uniCoords$PC2 - centroid[2]) ^ 2) } else {center <- NA} SppClims$nLocsCore[i] <- length(which(!is.na(MeanDists))) if (any(!is.na(MeanDists))) { if (nrow(uniCoords) >= 3) { SppClims$RhoCore[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate SppClims$RhoPCore[i] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } } SppClims$PC1.CMD <- SppSums$PC1.CMD SppClims$PC2.CMD <- SppSums$PC2.CMD SppClims$PC1.GM <- SppSums$PC1.GM SppClims$PC2.GM <- SppSums$PC2.GM SppClims$PC1.PC <- SppSums$PC1.PC SppClims$PC2.PC <- SppSums$PC2.PC SppClims$PC1max <- SppSums$PC1max SppClims$PC1min <- SppSums$PC1min SppClims$PC2max <- SppSums$PC2max SppClims$PC2min <- SppSums$PC2min ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Randomizations for geograhpic relationships ## ############################################################################################################################## ## Change year minimum here RData <- Data[which(Data$YearCount >= 2 & Data$YearsSurveyed >= 5), ] ######################################################################## # ## Randomizations Distance to Edge D2ERandG <- as.data.frame(matrix(NA, nrow = 1, ncol = 3)) names(D2ERandG) <- c('SppID', 'RhoCont', 'RhoPCont') for (i in 1:nrow(D2ESlopesG)) { current <- RData[which(RData$SppID == D2ESlopesG$SppID[i] & !is.na(RData$DistEdge)), ] if (nrow(current) > 0) { CurRands <- as.data.frame(matrix(NA, nrow = 100, ncol = 3)) names(CurRands) <- c('SppID', 'RhoCont', 'RhoPCont') CurRands$SppID <- D2ESlopesG$SppID[i] SpP <- SpatialPoints(current[, c('long', 'lat')], proj4string = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) SpP <- spTransform(SpP, CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) # Aggregate at cell level PRast <- rasterize(SpP, rbase, (current$Abun), fun = mean) DistRast <- rasterize(SpP, rbase, current$DistEdge, fun = mean) uniCoords <- xyFromCell(PRast, which(!is.na(values(PRast)) & !is.na(values(DistRast)))) ## Repeat process from gridded calculations 100 times for (k in 1:100) { # Shuffle abundance values current$Abun <- sample(current$Abun) # Aggregate at cell level PRast <- rasterize(SpP, rbase, (current$Abun), fun = mean) DistRast <- rasterize(SpP, rbase, current$DistEdge, fun = mean) uniCoords <- xyFromCell(PRast, which(!is.na(values(PRast)) & !is.na(values(DistRast)))) MeanBuns <- values(PRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast)))] MeanDists <- values(DistRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast)))] if (nrow(uniCoords) >= 3) { CurRands$RhoCont[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate CurRands$RhoPCont[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } D2ERandG <- rbind(D2ERandG, CurRands) } print(paste0(i, ' species plot-mapped')) } D2ERandG <- D2ERandG[2:nrow(D2ERandG), ] ### Random for centers GData <- Data[which(Data$YearCount >= 2 & Data$YearsSurveyed >= 5), ] GRandsG <- as.data.frame(matrix(NA, nrow = 1, ncol = 10)) names(GRandsG) <- c('SppId', 'Common.name', 'RhoPPC', 'RhoPC', 'RhoPGM', 'RhoGM', 'RhoPCMD', 'RhoCMD', 'RhoPAbun', 'RhoAbun') NAids <- c(132, 118, 61, 4, 3, 20, 182, 101, 153, 353, 411) for (i in 1:nrow(GSlopesG)) { if (GSlopesG$SppId[i] %in% NAids == F) { CurRands <- as.data.frame(matrix(NA, nrow = 100, ncol = 10)) names(CurRands) <- c('SppId', 'Common.name', 'RhoPPC', 'RhoPC', 'RhoPGM', 'RhoGM', 'RhoPCMD', 'RhoCMD', 'RhoPAbun', 'RhoAbun') CurRands$SppId <- GSlopesG$SppId[i] CurRands$Common.name <- GSlopesG$Common.name[i] current <- GData[which(GData$SppID == CurRands$SppId[1]), ] if (nrow(current) > 1) { SpP <- SpatialPoints(current[, c('long', 'lat')], proj4string = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) SpP <- spTransform(SpP, CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) PRast <- rasterize(SpP, rbase, (current$Abun), fun = mean) DistRast1 <- rasterize(SpP, rbase, current$DistPC, fun = mean) MeanBuns <- values(PRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast1)))] # Repeat 100 times for (k in 1:100) { # Shuffle abundances current$Abun <- sample(current$Abun) # Aggregate at cell level PRast <- rasterize(SpP, rbase, (current$Abun), fun = mean) DistRast1 <- rasterize(SpP, rbase, current$DistPC, fun = mean) DistRast2 <- rasterize(SpP, rbase, current$DistGM, fun = mean) DistRast3 <- rasterize(SpP, rbase, current$DistCMD, fun = mean) uniCoords <- xyFromCell(PRast, which(!is.na(values(PRast)) & !is.na(values(DistRast1)))) uniCoordsLL <- SpatialPoints(uniCoords[, 1:2], proj4string = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9')) uniCoordsLL <- spTransform(uniCoordsLL, CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) MeanBuns <- values(PRast)[which(!is.na(values(PRast)) & !is.na(values(DistRast1)))] MeanPC <- values(DistRast1)[which(!is.na(values(PRast)) & !is.na(values(DistRast1)))] MeanGM <- values(DistRast2)[which(!is.na(values(PRast)) & !is.na(values(DistRast2)))] MeanCMD <- values(DistRast3)[which(!is.na(values(PRast)) & !is.na(values(DistRast3)))] Top10 <- MeanBuns[which(MeanBuns >= quantile(MeanBuns, .9))] TopCoords <- uniCoords[which(MeanBuns >= quantile(MeanBuns, .9)), ] if (length(Top10) >= 1) { if (length(Top10) >= 2) { centroid <- gCentroid(SpatialPoints(TopCoords[, 1:2], proj4string = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'))) } else { centroid <- gCentroid(SpatialPoints(matrix(TopCoords[1:2], nrow = 1, ncol = 2), proj4string = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'))) } centerLL <- spTransform(centroid, CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) } if (nrow(uniCoords) >= 3 & !is.na(current$DistPC[1])) { CurRands$RhoPPC[k] <- cor.test(log10(MeanBuns), MeanPC, method = 'spearman')$p.value CurRands$RhoPC[k] <- cor.test(log10(MeanBuns), MeanPC, method = 'spearman')$estimate } if (nrow(uniCoords) >= 3 & !is.na(current$DistGM[1])) { CurRands$RhoPGM[k] <- cor.test(log10(MeanBuns), MeanGM, method = 'spearman')$p.value CurRands$RhoGM[k] <- cor.test(log10(MeanBuns), MeanGM, method = 'spearman')$estimate } if (nrow(uniCoords) >= 3 & !is.na(current$DistCMD[1])) { CurRands$RhoPCMD[k] <- cor.test(log10(MeanBuns), MeanCMD, method = 'spearman')$p.value CurRands$RhoCMD[k] <- cor.test(log10(MeanBuns), MeanCMD, method = 'spearman')$estimate } if (nrow(uniCoords) >= 3 & !is.na(centerLL[1])) { dists <- rdist.earth(centerLL@coords, current[, 4:3]) DistRast4 <- rasterize(SpP, rbase, dists[1, ], fun = mean) MeanCore <- values(DistRast4)[which(!is.na(values(PRast)) & !is.na(values(DistRast4)))] CurRands$RhoPAbun[k] <- cor.test(log10(MeanBuns), MeanCore, method = 'spearman')$p.value CurRands$RhoAbun[k] <- cor.test(log10(MeanBuns), MeanCore, method = 'spearman')$estimate } } GRandsG <- rbind(GRandsG, CurRands) } } print(paste0(i, ' species plot-mapped')) } GRandsG <- GRandsG[2:nrow(GRandsG), ] ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Randomizations for climate relationships ## ############################################################################################################################## ## vv change for min counts RClimData <- ClimData[which(ClimData$YearCount >= 2 & ClimData$YearsSurveyed >= 5), ] RandClimsG<- as.data.frame(matrix(NA, nrow = 1, ncol = 12)) names(RandClimsG) <- c('SppID', 'Common.name', 'RhoE', 'RhoPE', 'RhoPC', 'RhoPPC', 'RhoGM', 'RhoPGM', 'RhoCMD', 'RhoPCMD', 'RhoCore', 'RhoPCore') for (i in 1:nrow(SppClimsG)) { current <- RClimData[which(RClimData$SppID == SppClimsG$SppID[i] & !is.na(RClimData$dEdge)), ] if (nrow(current) > 1) { CurRands <- as.data.frame(matrix(NA, nrow = 100, ncol = 12)) names(CurRands) <- c('SppID', 'Common.name', 'RhoE', 'RhoPE', 'RhoPC', 'RhoPPC', 'RhoGM', 'RhoPGM', 'RhoCMD', 'RhoPCMD', 'RhoCore', 'RhoPCore') CurRands$SppID <- SppClimsG$SppID[i] CurRands$Common.name <- SppClimsG$Common.name[i] # grid (round) lat lon coordinates current$PC1 <- round_any(current$PC1, .2) current$PC2 <- round_any(current$PC2, .2) uniCoords <- unique(current[, c('PC1', 'PC2')]) # Repeat randomization 100 times for (k in 1:100) { # Shuffle abundance values current$Abun <- sample(current$Abun) MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dEdge[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } if (nrow(uniCoords) >= 3) { CurRands$RhoE[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate CurRands$RhoPE[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } ## Centers current <- RClimData[which(RClimData$SppID == SppClimsG$SppID[i] & !is.na(RClimData$dCMD)), ] # grid (round) lat lon coordinates current$PC1 <- round_any(current$PC1, .2) current$PC2 <- round_any(current$PC2, .2) uniCoords <- unique(current[, c('PC1', 'PC2')]) for (k in 1:100) { current$Abun <- sample(current$Abun) # PC MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dPC[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } if (nrow(uniCoords) >= 3) { CurRands$RhoPC[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate CurRands$RhoPPC[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } # GM MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dGM[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } if (nrow(uniCoords) >= 3) { CurRands$RhoGM[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate CurRands$RhoPGM[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } # CMD MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) MeanDists[j] <- mean(current$dCMD[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } if (nrow(uniCoords) >= 3) { CurRands$RhoCMD[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate CurRands$RhoPCMD[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } ## Core current <- RClimData[which(RClimData$SppID == SppClimsG$SppID[i]), ] # grid (round) lat lon coordinates current$PC1 <- round_any(current$PC1, .2) current$PC2 <- round_any(current$PC2, .2) uniCoords <- unique(current[, c('PC1', 'PC2')]) CurRands$nLocsCore <- nrow(uniCoords) for (k in 1:100) { current$Abun <- sample(current$Abun) MeanBuns <- NA MeanDists <- NA for (j in 1:nrow(uniCoords)) { MeanBuns[j] <- mean(current$Abun[which(current$PC1 == uniCoords$PC1[j] & current$PC2 == uniCoords$PC2[j])], na.rm = T) } Top10 <- MeanBuns[which(MeanBuns >= quantile(MeanBuns, .9, na.rm = T))] TopCoords <- uniCoords[which(MeanBuns >= quantile(MeanBuns, .9, na.rm = T)), ] if (length(Top10) >= 1) { centroid <- c(mean(TopCoords$PC1), mean(TopCoords$PC2)) CurRands[k, c('CorePC1', 'CorePC2')] <- centroid MeanDists <- sqrt((uniCoords$PC1 - centroid[1]) ^ 2 + (uniCoords$PC2 - centroid[2]) ^ 2) } else {center <- NA} if (any(!is.na(MeanDists))) { if (nrow(uniCoords) >= 3) { CurRands$RhoCore[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$estimate CurRands$RhoPCore[k] <- cor.test(log10(MeanBuns), MeanDists, method = 'spearman', use = 'complete')$p.value } } } RandClimsG <- rbind(RandClimsG, CurRands) } } RandClimsG <- RandClimsG[2:nrow(RandClimsG), ] ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Calculate coverage of species ranges by BBS surveys ## ############################################################################################################################## ## Get species that were in our analysis AllSpp <- unique(GSlopesG$SppId) ## Set up to capture data Coverage <- as.data.frame(matrix(NA, nrow = length(AllSpp), ncol = 14)) names(Coverage) <- c('SppID', 'MapName', 'Common.name', 'Order', 'Family', 'Area1', 'Cover1', 'PercentCov1', 'Area25', 'Cover25', 'PercentCov25', 'Area50', 'Cover50', 'PercentCov50') Coverage$SppID <- AllSpp ## Fill in other taxonomic info SppList$Order <- as.character(SppList$Order) SppList$Family <- as.character(SppList$Family) SppList$Common.Name <- as.character(SppList$Common.Name) SppList$scinamemap <- as.character(SppList$scinamemap) for (i in 1:nrow(Coverage)) { current <- SppList[which(SppList$SppID == Coverage$SppID[i]), ] if (nrow(current) > 0) { Coverage$MapName[i] <- current$scinamemap[1] Coverage$Order[i] <- current$Order[1] Coverage$Family[i] <- current$Family[1] Coverage$Common.name[i] <- current$Common.Name[1] } } # Set up for different minimum number of observations Coverage1 <- Coverage5 <- Coverage ## Create reference raster EqArea <- CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9') World <- getMap() World <- spTransform(World, EqArea) rbase <- raster(xmn = -165, xmx = -30, ymn = -60, ymx = 80, resolution = 1) # Three reference rasters correspond to three different resolutions for determining range coverage rbase1 <- projectRaster(rbase, crs = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'), res = c(100000,100000)) rbase25 <- projectRaster(rbase, crs = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'), res = c(250000,250000)) rbase50 <- projectRaster(rbase, crs = CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9'), res = c(500000,500000)) # gWorld <- gBuffer(World, byid = T, width = 50000) WorldR <- rasterize(World, rbase1, 'REGION') # gWorld <- gBuffer(World, byid = T, width = 125000) WorldR25 <- rasterize(World, rbase25, 'REGION') # gWorld <- gBuffer(World, byid = T, width = 250000) WorldR50 <- rasterize(World, rbase50, 'REGION') ## Loop to calculate coverage for (i in 1:nrow(Coverage)) { ## Check if map that fits scientific name exists if(file.exists(paste0("/Users/trvr/Box Sync/Bird Data/Bird Ranges/All/", Coverage$MapName[i], ".shp"))){ ## Open map map <- readOGR('/Users/trvr/Box Sync/Bird Data/Bird Ranges/All', layer = Coverage$MapName[i], verbose = F) ## Remove bits of range that aren't needed map <- map[map$PRESENCE != 3 & map$PRESENCE != 4 & map$PRESENCE != 5 & map$PRESENCE != 6 & map$SEASONAL != 3 & map$SEASONAL != 5 & map$SEASONAL != 4, ] #map <- unionSpatialPolygons(map, IDs = map$SISID) ## Reproject map to equal area EqArea <- CRS('+proj=laea +lat_0=30.2 +lon_0=-106.9') map_un <- map map <- spTransform(map, EqArea) MapR <- rasterize(map, rbase1) MapR25 <- rasterize(map, rbase25) MapR50 <- rasterize(map, rbase50) # Remove Non north and south american regions values(MapR)[-which(values(WorldR) %in% 6:7)] <- NA values(MapR25)[-which(values(WorldR25) %in% 6:7)] <- NA values(MapR50)[-which(values(WorldR50) %in% 6:7)] <- NA ## Calculates area (in # of gridcells) for each resolution Coverage$Area1[i] <- Coverage1$Area1[i] <- Coverage5$Area1[i] <- length(which(!is.na(values(MapR)))) Coverage$Area25[i] <- Coverage1$Area25[i] <- Coverage5$Area25[i] <- length(which(!is.na(values(MapR25)))) Coverage$Area50[i] <- Coverage1$Area50[i] <- Coverage5$Area50[i] <- length(which(!is.na(values(MapR50)))) ############### ## Pull out location data for species (min 2) current <- Data[which(Data$SppID == Coverage$SppID[i] & Data$YearCount >=2 & Data$YearsSurveyed >= 5), ] current$SCINAME <- 'Spp' if (nrow(current) > 0) { ## Create spatial points of locations, need to provide lon and lat pts <- SpatialPointsDataFrame(current[,c('long','lat')], proj4string = CRS(proj4string(map_un)), data = current) ## Reproject to match raster pts <- spTransform(pts, EqArea) ## Rasterize at each resolution PtsR <- rasterize(pts, MapR, 'SppID') PtsR25 <- rasterize(pts, MapR25, 'SppID') PtsR50 <- rasterize(pts, MapR50, 'SppID') ## Calculate coverage of points in absolute and %s Coverage$Cover1[i] <- length(which(!is.na(values(PtsR)) & !is.na(values(MapR)))) Coverage$PercentCov1[i] <- Coverage$Cover1[i]/Coverage$Area1[i] Coverage$Cover25[i] <- length(which(!is.na(values(PtsR25)) & !is.na(values(MapR25)))) Coverage$PercentCov25[i] <- Coverage$Cover25[i]/Coverage$Area25[i] Coverage$Cover50[i] <- length(which(!is.na(values(PtsR50)) & !is.na(values(MapR50)))) Coverage$PercentCov50[i] <- Coverage$Cover50[i]/Coverage$Area50[i] } ############### ## Pull out location data for species (min 1) current <- Data[which(Data$SppID == Coverage$SppID[i] & Data$YearsSurveyed >= 5), ] current$SCINAME <- 'Spp' if (nrow(current) > 0) { ## Create spatial points of locations, need to provide lon and lat pts <- SpatialPointsDataFrame(current[,c('long','lat')], proj4string = CRS(proj4string(map_un)), data = current) ## Reproject to match raster pts <- spTransform(pts, EqArea) ## Rasterize at each resolution PtsR <- rasterize(pts, MapR, 'SppID') PtsR25 <- rasterize(pts, MapR25, 'SppID') PtsR50 <- rasterize(pts, MapR50, 'SppID') ## Calculate coverage of points in absolute and %s Coverage1$Cover1[i] <- length(which(!is.na(values(PtsR)) & !is.na(values(MapR)))) Coverage1$PercentCov1[i] <- Coverage1$Cover1[i]/Coverage1$Area1[i] Coverage1$Cover25[i] <- length(which(!is.na(values(PtsR25)) & !is.na(values(MapR25)))) Coverage1$PercentCov25[i] <- Coverage1$Cover25[i]/Coverage1$Area25[i] Coverage1$Cover50[i] <- length(which(!is.na(values(PtsR50)) & !is.na(values(MapR50)))) Coverage1$PercentCov50[i] <- Coverage1$Cover50[i]/Coverage1$Area50[i] } ############### ## Pull out location data for species (min 5) current <- Data[which(Data$SppID == Coverage$SppID[i] & Data$YearCount >=5 & Data$YearsSurveyed >= 5), ] current$SCINAME <- 'Spp' if (nrow(current) > 0) { ## Create spatial points of locations, need to provide lon and lat pts <- SpatialPointsDataFrame(current[,c('long','lat')], proj4string = CRS(proj4string(map_un)), data = current) ## Reproject to match raster pts <- spTransform(pts, EqArea) ## Rasterize at each resolution PtsR <- rasterize(pts, MapR, 'SppID') PtsR25 <- rasterize(pts, MapR25, 'SppID') PtsR50 <- rasterize(pts, MapR50, 'SppID') ## Calculate coverage of points in absolute and %s Coverage5$Cover1[i] <- length(which(!is.na(values(PtsR)) & !is.na(values(MapR)))) Coverage5$PercentCov1[i] <- Coverage5$Cover1[i]/Coverage5$Area1[i] Coverage5$Cover25[i] <- length(which(!is.na(values(PtsR25)) & !is.na(values(MapR25)))) Coverage5$PercentCov25[i] <- Coverage5$Cover25[i]/Coverage5$Area25[i] Coverage5$Cover50[i] <- length(which(!is.na(values(PtsR50)) & !is.na(values(MapR50)))) Coverage5$PercentCov50[i] <- Coverage5$Cover50[i]/Coverage5$Area50[i] } } rm(map, gmap, map_un) gc() } ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Asess support for abundant-center and abundant-core predictions ##Geographic space ############################################################################################################################## ## In this section, ## Coverage threshold can be set to different values, or the datatables ## can be changed to check different criteria (including the non-gridded or Random tables) # ResTab Stores the proportion of species: # Column 1: With relationships in the predicted direction and > 10% of variation in abundance explained by range position # Column 2: Significant relationships in the predicted direction # Column 3: Non-significant relationships # Column 4: Significant relationships opposite to predictions # Column 5: With relationships opposite to predictions and > 10% of variation in abundance explained by range position # Column 6: the number of species in a given sample # Rows correspond to: # Geographic range edge - all birds # Geographic center - all birds # Goegraphic Abundance Core - all birds # Goegraphic range edge - Passerines # Goegraphic center - Passerines # Goegraphic Abundance Core - Passerines ResTab <- matrix(NA, nrow = 6, ncol = 6) # D2E GRID corz <- D2ESlopesG$RhoCont[which(Coverage$PercentCov1 >= .5)] pz <- D2ESlopesG$RhoPCont[which(Coverage$PercentCov1 >= .5)] ResTab[1,1] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,2] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[1,4] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,5] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,6] <- length(which(!is.na(corz))) # Center GRID corz <- GSlopesG$RhoCMD[which(Coverage$PercentCov1 >= .5)] pz <- GSlopesG$RhoPCMD[which(Coverage$PercentCov1 >= .5)] ResTab[2,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[2,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,6] <- length(which(!is.na(corz))) # Core GRID corz <- GSlopesG$RhoAbun[which(Coverage$PercentCov1 >= .5)] pz <- GSlopesG$RhoPAbun[which(Coverage$PercentCov1 >= .5)] ResTab[3,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[3,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,6] <- length(which(!is.na(corz))) # D2E Passerines GRID corz <- D2ESlopesG$RhoCont[which(Coverage$PercentCov1 >= .5 & GSlopesG$Order == 'Passeriformes')] pz <- D2ESlopesG$RhoPCont[which(Coverage$PercentCov1 >= .5 & GSlopesG$Order == 'Passeriformes')] ResTab[4,1] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,2] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[4,4] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,5] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,6] <- length(which(!is.na(corz))) # Center Passerines GRID corz <- GSlopesG$RhoCMD[which(Coverage$PercentCov1 >= .5 & GSlopesG$Order == 'Passeriformes')] pz <- GSlopesG$RhoPCMD[which(Coverage$PercentCov1 >= .5 & GSlopesG$Order == 'Passeriformes')] ResTab[5,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[5,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,6] <- length(which(!is.na(corz))) # Core Passerines GRID corz <- GSlopesG$RhoAbun[which(Coverage$PercentCov1 >= .5 & GSlopesG$Order == 'Passeriformes')] pz <- GSlopesG$RhoPAbun[which(Coverage$PercentCov1 >= .5 & GSlopesG$Order == 'Passeriformes')] ResTab[6,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[6,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,6] <- length(which(!is.na(corz))) ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Asess support for abundant-center and abundant-core predictions ## Climate space ############################################################################################################################## ## In this section, ## Coverage threshold can be set to different values, or the datatables ## can be changed to check different criteria (including the non-gridded or Random tables) # ResTab Stores the proportion of species: # Column 1: With relationships in the predicted direction and > 10% of variation in abundance explained by range position # Column 2: Significant relationships in the predicted direction # Column 3: Non-significant relationships # Column 4: Significant relationships opposite to predictions # Column 5: With relationships opposite to predictions and > 10% of variation in abundance explained by range position # Column 6: the number of species in a given sample # Rows correspond to: # Climate range edge - all birds # Climate center - alla birds # Climate Abundance Core - all birds # Climate range edge - Passerines # Climate center - Passerines # Climate Abundance Core - Passerines Covs <- Coverage$SppID[which(Coverage$PercentCov1 >= 0.5)] Pass <- SppSums$SppID[which(SppSums$Order == 'Passeriformes')] ResTab <- matrix(NA, nrow = 6, ncol = 6) # D2E GRID corz <- SppClimsG$RhoE[which(SppClimsG$SppID %in% Covs)] pz <- SppClimsG$RhoPE[which(SppClimsG$SppID %in% Covs)] ResTab[1,1] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,2] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[1,4] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,5] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[1,6] <- length(which(!is.na(corz))) # Center GRID corz <- SppClimsG$RhoCMD[which(SppClimsG$SppID %in% Covs)] pz <- SppClimsG$RhoPCMD[which(SppClimsG$SppID %in% Covs)] ResTab[2,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[2,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[2,6] <- length(which(!is.na(corz))) # Core GRID corz <- SppClimsG$RhoCore[which(SppClimsG$SppID %in% Covs)] pz <- SppClimsG$RhoPCore[which(SppClimsG$SppID %in% Covs)] ResTab[3,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[3,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[3,6] <- length(which(!is.na(corz))) # D2E Passerines GRID corz <- SppClimsG$RhoE[which(SppClimsG$SppID %in% Covs & SppClimsG$SppID %in% Pass)] pz <- SppClimsG$RhoPE[which(SppClimsG$SppID %in% Covs & SppClimsG$SppID %in% Pass)] ResTab[4,1] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,2] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[4,4] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,5] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[4,6] <- length(which(!is.na(corz))) # Center Passerines GRID corz <- SppClimsG$RhoCMD[which(SppClimsG$SppID %in% Covs & SppClimsG$SppID %in% Pass)] pz <- SppClimsG$RhoPCMD[which(SppClimsG$SppID %in% Covs & SppClimsG$SppID %in% Pass)] ResTab[5,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[5,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[5,6] <- length(which(!is.na(corz))) # Core Passerines GRID corz <- SppClimsG$RhoCore[which(SppClimsG$SppID %in% Covs & SppClimsG$SppID %in% Pass)] pz <- SppClimsG$RhoPCore[which(SppClimsG$SppID %in% Covs & SppClimsG$SppID %in% Pass)] ResTab[6,1] <- length(which(corz < -1*sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,2] <- length(which(corz < 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,3] <- length(which(pz > 0.05))/length(which(!is.na(corz))) ResTab[6,4] <- length(which(corz > 0 & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,5] <- length(which(corz > sqrt(.1) & pz < 0.05))/length(which(!is.na(corz))) ResTab[6,6] <- length(which(!is.na(corz))) ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Prepare data for the geography of deviations analyses ############################################################################################################################## # Get the latitudinal and longitudinal extent fo species ranges GeoExtents <- data.frame(SppID = GSlopesG$SppId) GeoExtents$Area <- GeoExtents$LatMax <- GeoExtents$LatMin <- GeoExtents$LonMax <- GeoExtents$LonMin <- NA for (i in 1:nrow(GSlopesG)) { sciname <- SppList$scinamemap[which(SppList$SppID == GSlopesG$SppId[i])][1] if(file.exists(paste0("Range Maps/", sciname, ".shp"))) { mapdata <- readOGR("Range Maps/", layer = sciname, verbose = F) map <-mapdata[mapdata$PRESENCE != 3 & mapdata$PRESENCE != 4 & mapdata$PRESENCE != 5 & mapdata$PRESENCE != 6 & mapdata$SEASONAL != 3 & mapdata$SEASONAL != 5 & mapdata$SEASONAL != 4, ] GeoExtents[i, 2:5] <- extent(map)[1:4] GeoExtents$Area[i] <- sum(area(map)) } print(paste0(i, ' species extented!')) } ################################################ ## Create data table for comparing geograhpic to climate stats UniSpp <- unique(ClimData$SppID) EvE <- as.data.frame(matrix(NA, nrow = length(UniSpp), ncol = 42)) names(EvE) <- c('SppID', 'Common.name', 'rho', 'rhop', 'geoRho', 'geoPRho', 'geoRhoCent', 'geoRhoPCent', 'geoRhoCore', 'geoRhoPCore', 'Reg', 'CoreLon', 'CoreLat', 'CoreD2Cont', 'climRho', 'climPRho', 'climRhoCent', 'climRhoPCent', 'climRhoCore', 'climRhoPCore', 'PC1.CMD', 'PC2.CMD', 'PC1.GM', 'PC2.GM', 'PC1.PC', 'PC2.PC', 'PC1min', 'PC1max', 'PC2min', 'PC2max', 'CorePC1', 'CorePC2', 'CentLon', 'CentLat', 'LonMin', 'LonMax', 'LatMin', 'LatMax', 'MeanPC1', 'SDPC1', 'MeanPC2', 'SDPC2') for (i in 1:length(UniSpp)) { EvE$SppID[i] <- UniSpp[i] EvE$Common.name[i] <- as.character(SppList$Common.Name[which(SppList$SppID == UniSpp[i])][1]) # Extract species x route data in climate and geo space # assure that they have the same routes in the same order curClim <- ClimData[which(ClimData$SppID == UniSpp[i] & ClimData$YearCount >= 1), ] curGeo <- Data[which(Data$SppID == UniSpp[i] & Data$RouteID %in% curClim$RouteID), ] curGeo <- curGeo[which(!is.na(curGeo$DistEdge)), ] curGeo <- curGeo[which(!is.na(curGeo$DistCMD)), ] curClim <- curClim[which(curClim$RouteID %in% curGeo$RouteID), ] curClim <- curClim[order(curClim$RouteID), ] curGeo <- curGeo[order(curGeo$RouteID), ] if (nrow(curClim) > 1 & nrow(curGeo) > 1) { # correspondence between climate and geographic range position curCor <- cor.test(curClim$dE, curGeo$DistEdge, method = 'spearman') EvE$rho[i] <- curCor$estimate EvE$rhop[i] <- curCor$p.value } # Fill in other previously calculated statistics # Geographic EvE$geoRho[i] <- D2ESlopesG$RhoCont[which(D2ESlopesG$SppID == UniSpp[i])] EvE$geoPRho[i] <- D2ESlopesG$RhoPCont[which(D2ESlopesG$SppID == UniSpp[i])] EvE$geoRhoCent[i] <- GSlopesG$RhoCMD[which(GSlopesG$SppId == UniSpp[i])] EvE$geoRhoPCent[i] <- GSlopesG$RhoPCMD[which(GSlopesG$SppId == UniSpp[i])] EvE$geoRhoCore[i] <- GSlopesG$RhoAbun[which(GSlopesG$SppId == UniSpp[i])] EvE$geoRhoPCore[i] <- GSlopesG$RhoPAbun[which(GSlopesG$SppId == UniSpp[i])] EvE$climRho[i] <- SppClimsG$RhoE[which(SppClimsG$SppID == UniSpp[i])] EvE$climPRho[i] <- SppClimsG$RhoPE[which(SppClimsG$SppID == UniSpp[i])] EvE$climRhoCent[i] <- SppClimsG$RhoCMD[which(SppClimsG$SppID == UniSpp[i])] EvE$climRhoPCent[i] <- SppClimsG$RhoPCMD[which(SppClimsG$SppID == UniSpp[i])] EvE$climRhoCore[i] <- SppClimsG$RhoCore[which(SppClimsG$SppID == UniSpp[i])] EvE$climRhoPCore[i] <- SppClimsG$RhoPCore[which(SppClimsG$SppID == UniSpp[i])] # Region of core EvE$Reg[i] <- D2ESlopesG$Reg[which(D2ESlopesG$SppID == UniSpp[i])] # Core position in space- lat, lon, distance to coast EvE$CoreLon[i] <- D2ESlopesG$CoreLon[which(D2ESlopesG$SppID == UniSpp[i])] EvE$CoreLat[i] <- D2ESlopesG$CoreLat[which(D2ESlopesG$SppID == UniSpp[i])] EvE$CoreD2Cont[i] <- D2ESlopesG$CorD2Cont[which(D2ESlopesG$SppID == UniSpp[i])] # Extent in scaled climate space EvE$PC1min[i] <- SppClimsG$PC1min[which(SppClimsG$SppID == UniSpp[i])] EvE$PC1max[i] <- SppClimsG$PC1max[which(SppClimsG$SppID == UniSpp[i])] EvE$PC2min[i] <- SppClimsG$PC2min[which(SppClimsG$SppID == UniSpp[i])] EvE$PC2max[i] <- SppClimsG$PC2max[which(SppClimsG$SppID == UniSpp[i])] # Center position in scaled climate space EvE$PC1.CMD[i] <- SppClimsG$PC1.CMD[which(SppClimsG$SppID == UniSpp[i])] EvE$PC2.CMD[i] <- SppClimsG$PC2.CMD[which(SppClimsG$SppID == UniSpp[i])] EvE$PC1.GM[i] <- SppClimsG$PC1.GM[which(SppClimsG$SppID == UniSpp[i])] EvE$PC2.GM[i] <- SppClimsG$PC2.GM[which(SppClimsG$SppID == UniSpp[i])] EvE$PC1.PC[i] <- SppClimsG$PC1.PC[which(SppClimsG$SppID == UniSpp[i])] EvE$PC2.PC[i] <- SppClimsG$PC2.PC[which(SppClimsG$SppID == UniSpp[i])] EvE$CorePC1[i] <- SppClimsG$CorePC1[which(SppClimsG$SppID == UniSpp[i])] EvE$CorePC2[i] <- SppClimsG$CorePC2[which(SppClimsG$SppID == UniSpp[i])] # cent lat and lon EvE$CentLong[i] <- Mids$x.CMD[which(Mids$SppID == UniSpp[i])] EvE$CentLat[i] <- Mids$y.CMD[which(Mids$SppID == UniSpp[i])] # lat and lon extent of range EvE$LonMin[i] <- GeoExtents$LonMin[which(GeoExtents$SppID == UniSpp[i])] EvE$LonMax[i] <- GeoExtents$LonMax[which(GeoExtents$SppID == UniSpp[i])] EvE$LatMin[i] <- GeoExtents$LatMin[which(GeoExtents$SppID == UniSpp[i])] EvE$LatMax[i] <- GeoExtents$LatMax[which(GeoExtents$SppID == UniSpp[i])] # Mean and standard deviation of climate from range # for converting from scaled to raw climate data EvE$MeanPC1[i] <- ClimStats$MeanPC1[which(ClimStats$SppID == UniSpp[i])] EvE$SDPC1[i] <- ClimStats$SDPC1[which(ClimStats$SppID == UniSpp[i])] EvE$MeanPC2[i] <- ClimStats$MeanPC2[which(ClimStats$SppID == UniSpp[i])] EvE$SDPC2[i] <- ClimStats$SDPC2[which(ClimStats$SppID == UniSpp[i])] } # Deviation of climate abundance core from climate center EvE$PC1dif <- EvE$CorePC1 - EvE$PC1.CMD EvE$PC2dif <- EvE$CorePC2 - EvE$PC2.CMD # Climate abundance core position converted to original temp and precip scale EvE$CorePC1 <- EvE$CorePC1*EvE$SDPC1 + EvE$MeanPC1 EvE$CorePC2 <- EvE$CorePC2*EvE$SDPC2 + EvE$MeanPC2 # Relative lat and lon deviations EvE$LatDif <- (EvE$CoreLat - EvE$CentLat)/(EvE$LatMax - EvE$LatMin) EvE$LonDif <- (EvE$CoreLon - EvE$CentLong)/(EvE$LonMax - EvE$LonMin) ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ############################################################################################################################## ## Test the geography of deviations ############################################################################################################################## # Restrict dataset based on coverage and taxa Covs <- Coverage$SppID[which(Coverage$PercentCov1 >= .5)] Pass <- SppList$SppID[which(SppList$Order == 'Passeriformes')] CurEvE <- EvE[which(EvE$SppID %in% Covs & EvE$SppID %in% Pass), ] # Proportion species with significant Climate-geography ranfe position relationship length(which(CurEvE$rho > 0 & CurEvE$rhop < 0.05))/nrow(CurEvE) # Test whether support for abundant-center is correlated with climate-geography correspondence # Geo space summary(lm(CurEvE$geoRho ~ CurEvE$rho)) # Test whether climate-geography correspondence varies among regions Anova(aov(CurEvE$rho ~ CurEvE$Reg)) # Test whether geographic abundant-center support varies among regions Anova(aov(CurEvE$geoRho ~ CurEvE$Reg)) # Test whether climate abundant-center support varies among regions Anova(aov(CurEvE$climRho ~ CurEvE$Reg)) # Test whether geographic abundant-center support increases with distance to coast summary(lm(CurEvE$geoRho ~ CurEvE$CoreD2Cont)) # Test whether climate abundant-center support increases with distance to coast summary(lm(CurEvE$climRho ~ CurEvE$CoreD2Cont)) # Test relationships between lat dif and latitude at core summary(lm(CurEvE$LatDif ~ CurEvE$CoreLat)) # Test relationships between lat dif and latitude at core summary(lm(CurEvE$LonDif ~ CurEvE$CoreLon)) # Test relationships between Temperature dif and temperature at core summary(lm(CurEvE$PC1dif ~ CurEvE$CorePC1)) # Test relationships between Precipitation dif and precipitation at core summary(lm(CurEvE$PC2dif ~ CurEvE$CorePC2))