This is the workflow in R of the published paper “Carranza, E., & Cardillo, M. (2024). The role of Magellanic penguins (Spheniscus magellanicus) on landscape dynamics and lithic taphonomy in La Pastosa Islet, North Patagonia, Argentina. The Holocene, 34(1), 87-97. https://doi.org/10.1177/09596836231200439”

##packages
library(geoR)
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-4 (built on 2024-02-14) is now loaded
## --------------------------------------------------------------
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 3.2-8
## Loading required package: spatstat.random
## spatstat.random 3.2-2
## Loading required package: spatstat.explore
## Loading required package: nlme
## spatstat.explore 3.2-6
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.2-10
## Loading required package: spatstat.linnet
## spatstat.linnet 3.1-4
## 
## spatstat 3.0-7 
## For an introduction to spatstat, type 'beginner'
library(sp)
library(maptools)
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: TRUE
## 
## Attaching package: 'maptools'
## The following object is masked from 'package:sp':
## 
##     sp2Mondrian
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggmap)
## Loading required package: ggplot2
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
library(ggplot2)
library(ggalluvial)

## base map
Pastosa_map<-get_googlemap(center = c(lon = -65.042, lat = -41.427),maptype = "satellite", zoom=16)
## ℹ <https://maps.googleapis.com/maps/api/staticmap?center=-41.427,-65.042&zoom=16&size=640x640&scale=2&maptype=satellite&key=xxx-0>
##plot map
ggmap(Pastosa_map)+
   xlab("Longitude") +  
  ylab("Latitude")

La pastosa Islet.Golfo San Matías. Rio Negro. Argentina

##data files   
cords<-read.table("cordPastosa.txt",T)# spatial coordinates
ocurr<-read.table("Ocurrence.txt",T)#ocurrence data
VC<-read.table("VC.txt",T)# % of vegetation cover
alt=read.table("Alture.txt",T)# Alture in m

##transform cualitative data into factors
ocurr$Artefacts=as.factor(ocurr$Artefacts)
ocurr$Nest=as.factor(ocurr$Nest)
ocurr$Ocurrence=as.factor(ocurr$Ocurrence)
#format and process spatial data
#read polygon
tkm <- getKMLcoordinates("PolPastosa.kml", ignoreAltitude=T)
## Warning: Code moved to https://github.com/rsbivand/spkml - seeking maintainer
## https://github.com/r-spatial/evolution/issues/6
##jitter redundant coordinates
jitcords <- jitterDupCoords(cords, max=0.00001)

##built a data frame 
esp1=data.frame(jitcords, ocurr)
head(esp1)
##get point patterns object with marks
complete_pp<- ppp(esp1$y,esp1$x, marks=esp1$Ocurrence, poly=tkm)
##plot densities for entire dataset
plot(density(split(complete_pp),sigma=0.0005), box=F)

##plot quadrant counts 
Quadrat <- quadratcount(split(complete_pp, nx= 6, ny=6))
plot(Quadrat, cex=0.5)

##Relative altitude for the area
al<- ppp(esp1$y,esp1$x, marks=alt, poly=tkm)
plot(Smooth(al, sigma=0.0005),box=F)
contour(Smooth(al, sigma=0.0005),add=T)

## vegetation cover map

cover1<- ppp(esp1$y,esp1$x, marks=VC, poly=tkm)
plot(Smooth(cover1, sigma=0.0005),box=F)
contour(Smooth(cover1, sigma=0.0005),add=T)

##filter by AP and PA to compare co-ocurrence pattern (see details in methods)
#AP (abscence of penguins, presence of artifacts) and
# PA (presence of penguins, abscence of artifacts)

esp2=esp1%>%
  filter(Ocurrence %in% c("AP", "PA"))

esp2$Ocurrence <- droplevels(esp2$Ocurrence)

##get point patterns object
APyPA_pp<- ppp(esp2$y,esp2$x, marks=esp2$Ocurrence, poly=tkm)


##plot densities for APyPA_pp(not shown in the paper)
plot(density(split(APyPA_pp),sigma=0.0005), box=F, main=NULL)

##plot quadrant counts form APyPA_pp(not shown in the paper)
Quadrat <- quadratcount(split(APyPA_pp, nx= 6, ny=6))
plot(Quadrat, cex=0.7)

##fit CSR model to multitype point pattern with two factor levels
# of AP (abscence of penguins, presence of artifacts) and
# PA (presence of penguins, abscence of artifacts)

fit1<-ppm(APyPA_pp~marks)
fit1
## Stationary multitype Poisson process
## Fitted to point pattern dataset 'APyPA_pp'
## 
## Possible marks: 'AP' and 'PA'
## 
## Log intensity:  ~marks
## 
## Intensities:
##  beta_AP  beta_PA 
##  2567837 81949408 
## 
##              Estimate      S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) 14.758574 0.1313064 14.501218 15.015930   *** 112.39795
## marksPA      3.463038 0.1333478  3.201681  3.724395   ***  25.96998
##fit K test over fitted model and estimate confidence intervals by means of 100
#simulations of complete spatial randomness

K<-envelope(fit1,verbose=F)
K
plot(K)

##segregation test between AP y PA
#suppressWarnings of likehood bandwith
SegTest1=suppressWarnings(segregation.test(APyPA_pp, nsim=100, verbose=FALSE))

SegTest1
## 
##  Monte Carlo test of spatial segregation of types
## 
## data:  APyPA_pp
## T = 57.551, p-value = 0.009901
##segregation test penguins  vs landcover
VC$VC=as.factor(VC$VC)##transform into factor
cover2<- ppp(esp1$y,esp1$x, marks=VC$VC, poly=tkm)
#suppressWarnings of likehood bandwith
SegTest2=suppressWarnings(segregation.test(cover2, nsim=100, verbose=FALSE))

SegTest2
## 
##  Monte Carlo test of spatial segregation of types
## 
## data:  cover2
## T = 927.62, p-value = 0.009901
#Artifact dataset with taphonomic variables
Artifact_Dat=read.table("Artifact_dat.txt",T)

Artifact_Dat
##plot artifact data with taphonomic variables
ggplot(data=Artifact_Dat,
       aes(axis1 = RawMat,  
           axis2 = Type,     
           axis3 = Frag,
           axis4 = Carbonation,
           axis5 = Alteration,
           y = Frec)) +
  geom_alluvium(aes(fill = RawMat)) +
  geom_stratum(width = 1/12, fill = "black", color = "grey") +
  geom_label(stat = "stratum", 
             aes(label = after_stat(stratum))) +
             scale_x_discrete(limits = c("Response"),
             expand = c(0.15, 0.05)) +
              theme_void()

Curated and published by Marcelo Cardillo and Eugenia Carranza (IMHICIHU-CONICET-UBA)