#install.packages("sp") #install.packages("rgeos") library(sp) library(rgeos) ################################################################################## # Calculates the unbiased density for each sampling units (e.g., coral colony) based on the coordinates of the sampling points # # Arguments: (input) # Coor - a data frame with the X and Y coordinates of the sampling points within a quadrat OR # a vector of locations on the X axis only for the linear point intercept method # Diameters - a vector containing the diameters of the sampled individuals # quadsegs - number of line segments to use to approximate a quarter circle (from gBUffer in rgeos). Can be reduced to increase speed or increased for increased accuracy. # # Value: a vector of the estimated density of each sampling unit. 1/ESA - effective sampled area (ESA). #Note - NA values are stripped from Coor, locations of NA within Diameters are returned as NA # Jonathan Belmaker # 18 March 2019 ################################################################################## Density <- function(Coor, Diameters , quadsegs = 10) { Coor = data.frame(Coor) DiametersNA = is.na(Diameters) # finds NAs Diameters[DiametersNA] = 1 # replaces NA with 1, just so the function can run DiametersZero = (Diameters == 0) # finds zeros Diameters[DiametersZero] = 1 # replaces zeros with 1, just so the function can run if (dim(Coor)[2] == 1 ) {#adding the second axis for the linear point intercept method Coor= cbind(Coor, rep(0, length(Coor))) } Coor = Coor[complete.cases(Coor),] #making sure to remove NA values colnames(Coor) = c("X","Y") Radii = Diameters/2 Sampling_Locations <- SpatialPointsDataFrame(Coor, Coor) #make spatial point data frame Buffers = lapply(1:length(Radii), function(n) gBuffer(Sampling_Locations, width=Radii[n], quadsegs = quadsegs ))#apply buffer ESA = sapply(1:length(Radii), function(n) gArea(Buffers[[n]])) #calculate area Density = (1/ESA) Density[DiametersNA] = NA # puts back NA values Density[DiametersZero] = NA # puts back NA values return(Density) } ###This is were the code for this specific data starts: The data contains observations (rows) and metadata. ### It is essential to have this columns: Individual diameter (m), Points intervals (m), # of points per transect. ###change the next two lines to fit your directors and file name setwd("C:/Users/Documents/Project/Results") Data= read.csv("rawdata_for_density.csv") diameterVector = (Data$`Individual diameter (m)` < 0.01) Data$`Individual diameter (m)`[diameterVector] = 0.01 # A 1 cm threshold. DensityVector = sapply(1:nrow(Data), function(n) Density(Coor = seq(0, by = Data$`Points intervals (m)`[n], length.out = Data$`# of points per transect`[n]), Diameters = Data$`Individual diameter (m)`[n]))#apply function for all rows Data$Density = DensityVector write.csv(Data, file = "calculated_density.csv" )