Penguin_Density_Calculations.R

######################################################################
#        		    Code from Francomano et al. (IBIS 2023) 		    	   #


######################################################################
#                   Penguin density calculations   	    	           #
######################################################################

# Created by:	Dante Francomano (dfrancomano@alumni.purdue.edu)

# Date: 20201206

######################################################################
######################################################################

# Load packages
library(openxlsx)
library(rgdal)
library(dplyr)
library(dismo)

# Read in data
recorderLocations<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Site_Coordinates.csv")
IDLEPMRecorderLocations<-recorderLocations[1:2,]
IDLEPPARecorderLocations<-recorderLocations[3:7,]
IMRecorderLocations<-recorderLocations[8:14,]
IDLEGridPoints<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Datos_Espaciales/Isla_de_los_Estados_Uncut_Grid_Points.csv")
IMPolygon<-readOGR("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Datos_Espaciales/Isla_Martillo_Polygon_20190203.kml")
IMNoBeachPolygon<-readOGR("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Datos_Espaciales/Isla_Martillo_Sin_Playa_Polygon_20190203.kml")
IMStakeLocations<-read.xlsx("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Datos_Espaciales/Martillo_Estacas_2015.xlsx")[,c(1,3,2)]
PMCountDataPaths<-list.files("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Pingüinos/Conteos_de_Magallánicos",full.names=TRUE)
rawLocalIDLEPMCountData<-lapply(PMCountDataPaths[5:6],read.xlsx)
rawLocalIMCountData<-lapply(PMCountDataPaths[7:13],read.xlsx)
rawGlobalIMCountData<-read.csv(PMCountDataPaths[1])
rawPPAPhotoCountData<-read.csv("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Pingüinos/Fotos_para_PPA/Conteos_de_PPA_para_Analisis.csv",na.strings="N/A")
splitPPAPhotoCountData<-split(rawPPAPhotoCountData,rawPPAPhotoCountData[,1])

# Create containers for results
PPAResults<-as.data.frame(matrix(NA,5,126))
colnames(PPAResults)<-c("Site","Penguins_per_Square_Meter_within_1_Meter",sapply(2:125,function(x) paste("Penguins_per_Square_Meter_within_",x,"_Meters",sep="")))
PPAResults[,1]<-IDLEPPARecorderLocations[,1]
PMResults<-as.data.frame(matrix(NA,9,126))
colnames(PMResults)<-c("Site","Active_Nests_per_Square_Meter_within_1_Meter",sapply(2:125,function(x) paste("Active_Nests_per_Square_Meter_within_",x,"_Meters",sep="")))
PMResults[,1]<-c(IDLEPMRecorderLocations[,1],IMRecorderLocations[,1])


######################################################################
# Clean data
######################################################################

# Remove approximate duplicate stakes and clearly unused stakes from table of locations and remove lettering to match global count naming
IMStakeLocationsNoDuplicates<-IMStakeLocations[-which(IMStakeLocations[,1]=="035M"|IMStakeLocations[,1]=="P34"|IMStakeLocations[,1]=="P53"|IMStakeLocations[,1]=="P71"|IMStakeLocations[,1]=="069M"),]
IMStakeLocationsNoDuplicatesNoWeird<-IMStakeLocationsNoDuplicates[-which(IMStakeLocationsNoDuplicates[,1]=="AGUA37"|IMStakeLocationsNoDuplicates[,1]=="AGUA56"|IMStakeLocationsNoDuplicates[,1]=="FAROMA"),]
IMStakeLocationsNoDuplicatesNoWeirdRenamed<-IMStakeLocationsNoDuplicatesNoWeird
IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1]<-gsub("P","",IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1])
IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1]<-gsub("M","",IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1])
IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1]<-gsub("UEVO","",IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1])

# Remove leading zeros to match global count naming and rename L00 to 300 to homogenize names
IMStakeLocationsNoDuplicatesNoWeirdRenamed[which(IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1]=="L00"),1]<-300
IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1]<-as.numeric(IMStakeLocationsNoDuplicatesNoWeirdRenamed[,1])
IMStakeLocationsNoDuplicatesNoWeirdRenamed<-arrange(IMStakeLocationsNoDuplicatesNoWeirdRenamed,Name)

# Rename L00 to 300 in global count naming and remove count data with an unknown stake
rawGlobalIMCountData[,3]<-as.character(rawGlobalIMCountData[,3])
rawGlobalIMCountData[which(rawGlobalIMCountData[,3]=="L00"),3]<-"300"
globalIMCountDataCleaned<-rawGlobalIMCountData[which(rawGlobalIMCountData[,3]!="???"),]
globalIMCountDataCleaned[,3]<-as.numeric(globalIMCountDataCleaned[,3])
globalIMCountDataCleaned<-arrange(globalIMCountDataCleaned,Estaca,Distancia)


######################################################################
# Determine, assign, and write PPA results
######################################################################

# Synthesize and arrange data
rawPPAPhotoCountData$Penguins_per_Square_Meter<-rawPPAPhotoCountData$Number_of_Penguins/(pi*20^2)
PPASitePhotoCountData<-rawPPAPhotoCountData[which(is.na(rawPPAPhotoCountData$Full_Grid_Point)),]
PPANonSitePhotoCountData<-rawPPAPhotoCountData[which(!is.na(rawPPAPhotoCountData$Full_Grid_Point)),]
PPASitePhotoCountData$Latitude<-IDLEPPARecorderLocations$Latitude
PPASitePhotoCountData$Longitude<-IDLEPPARecorderLocations$Longitude
PPANonSitePhotoCountData$Latitude<-NA
PPANonSitePhotoCountData$Longitude<-NA
for(i in 1:dim(PPANonSitePhotoCountData)[1]) PPANonSitePhotoCountData$Latitude[i]<-IDLEGridPoints$Latitude[which(IDLEGridPoints$Raw_Point_ID==PPANonSitePhotoCountData$Full_Grid_Point[i])]
for(i in 1:dim(PPANonSitePhotoCountData)[1]) PPANonSitePhotoCountData$Longitude[i]<-IDLEGridPoints$Longitude[which(IDLEGridPoints$Raw_Point_ID==PPANonSitePhotoCountData$Full_Grid_Point[i])]


# Create custom coordinate reference system, and project data to it
IDLECRS<-CRS("+proj=aea +lat_1=-54.88 +lat_2=-54.87 +lat_0=-54.895 +lon_0=-64.66 +x_0=2000 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
spatialIDLERecorderLocationsLL<-SpatialPointsDataFrame(IDLEPPARecorderLocations[,3:2],as.data.frame(IDLEPPARecorderLocations[,1]),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
spatialIDLERecorderLocationsAEA<-spTransform(spatialIDLERecorderLocationsLL,IDLECRS)
spatialSynthesizedIDLECountPointDataLL<-SpatialPointsDataFrame(PPANonSitePhotoCountData[,14:13],PPANonSitePhotoCountData[,c(1:3,11,12)],proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
spatialSynthesizedIDLECountPointDataAEA<-spTransform(spatialSynthesizedIDLECountPointDataLL,IDLECRS)
plot(spatialIDLERecorderLocationsAEA)
text(spatialSynthesizedIDLECountPointDataAEA,spatialSynthesizedIDLECountPointDataAEA@data$Full_Grid_Point,cex=0.2)

# Extend count areas with Voronoi polygons
voronoiIDLEPolygonData<-voronoi(spatialSynthesizedIDLECountPointDataAEA)
plot(voronoiIDLEPolygonData)

# Remove 20-m recorder radii from Voronoi polygons
IDLERecorder20MeterPolygonData<-buffer(spatialIDLERecorderLocationsAEA,20,dissolve=FALSE)
clippedVoronoiIDLEPolygonDataRecorderRadiiExtracted<-raster::erase(voronoiIDLEPolygonData,IDLERecorder20MeterPolygonData)

# Assign densities to recorder polygons, homogenize data structures, and bind layers
IDLERecorder20MeterPolygonData@data$Penguins_per_Square_Meter<-PPASitePhotoCountData$Penguins_per_Square_Meter
names(IDLERecorder20MeterPolygonData@data)<-c("Location","Penguins_per_Square_Meter")
clippedVoronoiIDLEPolygonDataRecorderRadiiExtracted@data<-data.frame(Location=clippedVoronoiIDLEPolygonDataRecorderRadiiExtracted@data$Full_Grid_Point,Penguins_per_Square_Meter=clippedVoronoiIDLEPolygonDataRecorderRadiiExtracted@data$Penguins_per_Square_Meter)
IDLEAllPolygonData<-bind(IDLERecorder20MeterPolygonData,clippedVoronoiIDLEPolygonDataRecorderRadiiExtracted)

# Find maximum density and create colors based on it, and plot to visually verify results
IDLEMaximum20MeterDensity<-max(IDLEAllPolygonData@data$Penguins_per_Square_Meter)
IDLEColors<-sapply(IDLEAllPolygonData@data$Penguins_per_Square_Meter, function(x) rgb(x/IDLEMaximum20MeterDensity,x/IDLEMaximum20MeterDensity,x/IDLEMaximum20MeterDensity))
plot(IDLEAllPolygonData,col=IDLEColors)

# Generate varying radius plots around recorders, determine penguin density within them, and write results
IDLEVaryingRadiusPlots<-sapply(1:125,function(x) buffer(spatialIDLERecorderLocationsAEA,x,dissolve=FALSE))
for(i in 1:125) names(IDLEVaryingRadiusPlots[[i]]@data)<-"Site"
IDLEvaryingRadiusPlotIntersections<-sapply(IDLEVaryingRadiusPlots,function(x) raster::intersect(IDLEAllPolygonData,x))
plot(IDLEvaryingRadiusPlotIntersections[[50]],lwd=0.2)
for(i in 1:125) IDLEvaryingRadiusPlotIntersections[[i]]@data$Area<-area(IDLEvaryingRadiusPlotIntersections[[i]])
for(i in 1:125) IDLEvaryingRadiusPlotIntersections[[i]]@data$Number_of_Penguins<-IDLEvaryingRadiusPlotIntersections[[i]]@data$Area*IDLEvaryingRadiusPlotIntersections[[i]]@data$Penguins_per_Square_Meter
for(i in 1:5) for(j in 1:125) PPAResults[i,(j+1)]<-sum(IDLEvaryingRadiusPlotIntersections[[j]]@data$Number_of_Penguins[which(IDLEvaryingRadiusPlotIntersections[[j]]@data$Site==PPAResults$Site[i])])/(pi*j^2)
write.csv(PPAResults,"/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Pingüinos/PPA_Densities_by_Radii.csv")


######################################################################
# Assign local PM results
######################################################################

for(i in 1:20) PMResults[1:2,(i+1)]<-sapply(rawLocalIDLEPMCountData,function(x) length(which(x[,seq(5,19,2)]<=i))/(pi*i^2))
for(i in 1:20) PMResults[3:9,(i+1)]<-sapply(rawLocalIMCountData,function(x) length(which(x[,seq(5,19,2)]<=i))/(pi*i^2))


######################################################################
# Determine, assign, and write global PM results
######################################################################

# Synthesize IM global count data, assign true zeros, remove zeros where there are no data
synthesizedGlobalIMCountData<-cbind(IMStakeLocationsNoDuplicatesNoWeirdRenamed,rep(NA,dim(IMStakeLocationsNoDuplicatesNoWeirdRenamed)[1]),rep(NA,dim(IMStakeLocationsNoDuplicatesNoWeirdRenamed)[1]))
colnames(synthesizedGlobalIMCountData)[4:5]<-c("Number_of_Nests","Nests_per_Square_Meter")
for(i in 1:dim(synthesizedGlobalIMCountData)[1]) synthesizedGlobalIMCountData$Number_of_Nests[i]<-ifelse((length(which(unique(globalIMCountDataCleaned$Estaca)==synthesizedGlobalIMCountData$Name[i]))>0),dim(globalIMCountDataCleaned[which(globalIMCountDataCleaned$Estaca==synthesizedGlobalIMCountData$Name[i]),])[1],NA)
trueZeroSites<-globalIMCountDataCleaned$Estaca[which(globalIMCountDataCleaned$Distancia=="0_Nidos_Activos")]
for(i in trueZeroSites) synthesizedGlobalIMCountData$Number_of_Nests[synthesizedGlobalIMCountData$Name==i]<-0
synthesizedGlobalIMCountData$Nests_per_Square_Meter<-synthesizedGlobalIMCountData$Number_of_Nests/(pi*20^2)
synthesizedGlobalIMCountDataNoNA<-synthesizedGlobalIMCountData[which(!is.na(synthesizedGlobalIMCountData$Number_of_Nests)),]

# Create custom coordinate reference system, and project data to it
IMCRS<-CRS("+proj=aea +lat_1=-54.906 +lat_2=-54.908 +lat_0=-54.907 +lon_0=-67.383 +x_0=1200 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
spatialIMRecorderLocationsLL<-SpatialPointsDataFrame(IMRecorderLocations[,3:2],as.data.frame(IMRecorderLocations[,1]),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
spatialIMRecorderLocationsAEA<-spTransform(spatialIMRecorderLocationsLL,IMCRS)
projectedIMPolygonAEA<-spTransform(IMPolygon,IMCRS)
projectedIMNoBeachPolygonAEA<-spTransform(IMNoBeachPolygon,IMCRS)
spatialSynthesizedGlobalIMCountPointDataLL<-SpatialPointsDataFrame(synthesizedGlobalIMCountDataNoNA[,3:2],synthesizedGlobalIMCountDataNoNA[,c(1,4,5)],proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
spatialSynthesizedGlobalIMCountPointDataAEA<-spTransform(spatialSynthesizedGlobalIMCountPointDataLL,IMCRS)
plot(projectedIMPolygonAEA)
plot(projectedIMNoBeachPolygonAEA,add=TRUE)
plot(spatialIMRecorderLocationsAEA,add=TRUE)
text(spatialSynthesizedGlobalIMCountPointDataAEA,spatialSynthesizedGlobalIMCountPointDataAEA@data$Name,cex=0.5)

# Extend count areas to the full island with Voronoi polygons
voronoiGlobalIMPolygonData<-voronoi(spatialSynthesizedGlobalIMCountPointDataAEA,ext=as.numeric(bbox(projectedIMPolygonAEA))[c(1,3,2,4)])
clippedVoronoiGlobalIMPolygonData<-raster::intersect(voronoiGlobalIMPolygonData,projectedIMNoBeachPolygonAEA)
plot(projectedIMPolygonAEA)
plot(clippedVoronoiGlobalIMPolygonData,add=TRUE)

# Remove 20-m recorder radii from Voronoi polygons
IMRecorder20MeterPolygonData<-buffer(spatialIMRecorderLocationsAEA,20,dissolve=FALSE)
clippedVoronoiGlobalIMPolygonDataRecorderRadiiExtracted<-raster::erase(clippedVoronoiGlobalIMPolygonData,IMRecorder20MeterPolygonData)

# Assign densities to recorder polygons, homogenize data structures, and bind layers
IMRecorder20MeterPolygonData@data$Nests_per_Square_Meter<-PMResults[3:9,21]
names(IMRecorder20MeterPolygonData@data)<-c("Location","Nests_per_Square_Meter")
clippedVoronoiGlobalIMPolygonDataRecorderRadiiExtracted@data<-data.frame(Location=clippedVoronoiGlobalIMPolygonDataRecorderRadiiExtracted@data$Name.1,Nests_per_Square_Meter=clippedVoronoiGlobalIMPolygonDataRecorderRadiiExtracted@data$Nests_per_Square_Meter)
IMAllPolygonData<-bind(IMRecorder20MeterPolygonData,clippedVoronoiGlobalIMPolygonDataRecorderRadiiExtracted)

# Find maximum density and create colors based on it, and plot to visually verify results
IMMaximum20MeterDensity<-max(IMAllPolygonData@data$Nests_per_Square_Meter)
IMColors<-sapply(IMAllPolygonData@data$Nests_per_Square_Meter, function(x) rgb(x/IMMaximum20MeterDensity,x/IMMaximum20MeterDensity,x/IMMaximum20MeterDensity))
plot(projectedIMPolygonAEA)
plot(IMAllPolygonData,add=TRUE,col=IMColors)

# Generate varying radius plots around recorders, determine penguin density within them, and write results
IMVaryingRadiusPlots<-sapply(21:125,function(x) buffer(spatialIMRecorderLocationsAEA,x,dissolve=FALSE))
for(i in 1:105) names(IMVaryingRadiusPlots[[i]]@data)<-"Site"
IMvaryingRadiusPlotIntersections<-sapply(IMVaryingRadiusPlots,function(x) raster::intersect(IMAllPolygonData,x))
plot(projectedIMPolygonAEA)
plot(IMvaryingRadiusPlotIntersections[[50]],add=TRUE,lwd=0.2)
for(i in 1:105) IMvaryingRadiusPlotIntersections[[i]]@data$Area<-area(IMvaryingRadiusPlotIntersections[[i]])
for(i in 1:105) IMvaryingRadiusPlotIntersections[[i]]@data$Number_of_Nests<-IMvaryingRadiusPlotIntersections[[i]]@data$Area*IMvaryingRadiusPlotIntersections[[i]]@data$Nests_per_Square_Meter
for(i in 1:7) for(j in 1:105) PMResults[(i+2),(j+21)]<-sum(IMvaryingRadiusPlotIntersections[[j]]@data[which(IMvaryingRadiusPlotIntersections[[j]]@data$Site==PMResults[(i+2),1]),]$Number_of_Nests)/(pi*(j+20)^2)
write.csv(PMResults,"/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Conteos_de_Pingüinos/PM_Densities_by_Radii.csv")

# Create a figure showing process
pdf("/Users/dantefrancomano/Documents/Purdue/Dissertation/Pingüinos/Figures/Density_Calculation_Process.pdf",width=3,height=8,pointsize=6)
par(mfcol=c(4,1),mar=c(2,0,0,0),oma=c(0,2,2,2))
plot(projectedIMPolygonAEA,bg=rgb(223/255,231/255,236/255),col=rgb(239/255,231/255,217/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(projectedIMNoBeachPolygonAEA,add=TRUE,col=rgb(225/255,232/255,219/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(spatialSynthesizedGlobalIMCountPointDataAEA,add=TRUE,pch=20,cex=0.25,col=rgb(62/255,62/255,62/255))
plot(buffer(spatialSynthesizedGlobalIMCountPointDataAEA,20,dissolve=FALSE),add=TRUE,border=rgb(62/255,62/255,62/255),lwd=0.5)
#scaleBar(IMCRS,pos="bottomleft",cex=2,legend="",seg.len=10,lwd=1,col="grey30",bty="n",outer=FALSE)
lines(c(70,570),rep(-450,2),col=rgb(153/255,153/255,153/255),lwd=1)
lines(rep(70,2),c(-450,-430),col=rgb(153/255,153/255,153/255),lwd=1)
lines(rep(320,2),c(-450,-430),col=rgb(153/255,153/255,153/255),lwd=1)
lines(rep(570,2),c(-450,-430),col=rgb(153/255,153/255,153/255),lwd=1)
text(c(70,320,640),rep(-350,3),c("0","250","500 m"),cex=3,col=rgb(153/255,153/255,153/255))
text(100,620,"A.",cex=4,col=rgb(62/255,62/255,62/255))
box(col=rgb(153/255,153/255,153/255),lwd=4)

plot(projectedIMPolygonAEA,bg=rgb(223/255,231/255,236/255),col=rgb(239/255,231/255,217/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(projectedIMNoBeachPolygonAEA,add=TRUE,col=rgb(225/255,232/255,219/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(spatialSynthesizedGlobalIMCountPointDataAEA,add=TRUE,pch=20,cex=0.25,col=rgb(62/255,62/255,62/255))
plot(clippedVoronoiGlobalIMPolygonData,add=TRUE,border=rgb(62/255,62/255,62/255),lwd=0.5)
text(100,620,"B.",cex=4,col=rgb(62/255,62/255,62/255))
box(col=rgb(153/255,153/255,153/255),lwd=4)

plot(projectedIMPolygonAEA,bg=rgb(223/255,231/255,236/255),col=rgb(239/255,231/255,217/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(projectedIMNoBeachPolygonAEA,add=TRUE,col=rgb(225/255,232/255,219/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(spatialIMRecorderLocationsAEA,add=TRUE,pch=20,cex=0.25,col=rgb(62/255,62/255,62/255))
plot(IMvaryingRadiusPlotIntersections[[30]],add=TRUE,border=rgb(62/255,62/255,62/255),lwd=0.5)
text(100,620,"C.",cex=4,col=rgb(62/255,62/255,62/255))
box(col=rgb(153/255,153/255,153/255),lwd=4)

plot(projectedIMPolygonAEA,bg=rgb(223/255,231/255,236/255),col=rgb(239/255,231/255,217/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(projectedIMNoBeachPolygonAEA,add=TRUE,col=rgb(225/255,232/255,219/255),border=rgb(62/255,62/255,62/255),lwd=0.5)
plot(spatialIMRecorderLocationsAEA,add=TRUE,pch=20,cex=0.25,col=rgb(62/255,62/255,62/255))
plot(IMvaryingRadiusPlotIntersections[[60]],add=TRUE,border=rgb(62/255,62/255,62/255),lwd=0.5)
text(100,620,"D.",cex=4,col=rgb(62/255,62/255,62/255))
box(col=rgb(153/255,153/255,153/255),lwd=4)
dev.off()