
# MAKE MAMMAL DIVERSITY DATASETS
# from cleaned genus occurrence files (Data_GenusOccurrence.txt)
# also uses the grid (Data_RegionsGridded.shp)

# continental (NA, EU) and regional (6 regions)
# with global time bins (stratigraphic stages) and continental time bins (land mammal ages)

# Line 30:  summarize datasets
# Line 65:  compile raw diversity data
# Line 170: apply sampling corrections:
# - make occurrence and abundance matrices for regions and continents
# - run sampling correction algorithms
# Line 530: add size of sampled area to the datasets, and add which grid cells each
#           locality is in to the GenusOccurrences dataset

# produces:
# Data_Regions and Data_Continents (files for diversity-productivity analyses,
#    included in upload as changed subsequent to this by Rscript_ClimateDataProcessing.r)

# modifies Data_GenusOccurrence.txt as included in upload (adds grid cell ID for each locality)


# Script by Susanne Fritz, last changed 5th Oct 2015

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


# Summarize the occurrence dataset

# set your working directory!
setwd("D://Data//Gammadiv//FilesForUpload")
# NOTE: the following code will overwrite some datasets included in the upload.
# To avoid this, change file names.

ds <- read.delim("Data_GenusOccurrence.txt", as.is=TRUE) # load the occurrence data
head(ds)
any(duplicated(ds)) # any data lines duplicated? should be FALSE.
ds$latlong <- paste(ds$LAT, ds$LONG, sep="_") # spatially unique location identifiers
unique(ds$Continent)
summary(as.factor(ds$Region))
luni <- function(x) (length(unique(x))) # function to count unique entries in a vector
length(unique(ds$Genus))
length(unique(ds$Family))
tapply(ds$Genus, ds$Continent, luni) # how many genera in each continent
tapply(ds$Genus, ds$SourceDataset, luni) # how many in each source dataset
tapply(ds$latlong, ds$Continent, luni) # how many locations are unique for each continent?
tapply(ds$latlong, ds$SourceDataset, luni)
tapply(ds$Genus, ds$Continent, length) # how many genera in each continent?
tapply(ds$Genus, ds$Order, luni) # how many genera in each order?
tapply(ds$latlong, ds$Order, luni) # how many locations 
# more data summary for Extended Data Table 1
unigltime <- unique(ds$GlobalTimeBins)
for (i in seq(along=unigltime)) {
  print(unigltime[i])
  dsi <- ds[which(ds$GlobalTimeBins==unigltime[i]),]
  print(tapply(dsi$Genus, dsi$Order, luni))
  print(tapply(dsi$latlong, dsi$Order, luni))
}


# raw gamma diversity without sampling correction


# genus dataset for the regions, with global and continental time bins

unireg <- unique(ds$Region) # the regions
unicotime <- unique(ds$ContinentalBins[!is.na(ds$ContinentalBins)]) # the land mammal ages / continental time bins
unicotime <- unicotime[c(1:18,25,26,19:23,27:30,24,31)] # put into chronological order (first Europe then NA)
ds$ContinentalBins <- as.factor(ds$ContinentalBins)
ds$GlobalTimeBins <- as.factor(ds$GlobalTimeBins)
for (i in seq(along=unireg)) { # loop through the regions
  # global time bins
  dsi <- ds[which(ds$Region==unireg[i]),] # dataset for region i
  numgen <- tapply(dsi$Genus, dsi$GlobalTimeBins, luni) # total number of genera within each global time bin in region i
  numloc <- tapply(dsi$latlong, dsi$GlobalTimeBins, luni) # total number of localities within each global time bin in region i
  # combine the data together to a new dataset
  newline <- data.frame(Region=rep(unireg[i], length(unigltime)), GlobalTimeBins=unigltime, ContinentalBins=rep(999, length(unigltime)), NumLocs=numloc[match(unigltime, names(numloc))], NumGenera=numgen[match(unigltime, names(numgen))])
  newline$GlobalTimeBins <- as.character(newline$GlobalTimeBins)
  # continental time bins
  dsi <- ds[which(ds$Region==unireg[i] & !is.na(ds$ContinentalBins)),] # dataset for region i
  numgen <- tapply(dsi$Genus, dsi$ContinentalBins, luni) # total number of genera within each continental time bin in region i
  numloc <- tapply(dsi$latlong, dsi$ContinentalBins, luni) # total number of localities within each continental time bin in region i
  # combine the data together to a new dataset
  newline2 <- data.frame(Region=rep(unireg[i], length(unicotime)), GlobalTimeBins=rep(999, length(unicotime)), ContinentalBins=unicotime, NumLocs=numloc[match(unicotime, names(numloc))], NumGenera=numgen[match(unicotime, names(numgen))])
  newline2$ContinentalBins <- as.character(newline2$ContinentalBins)
  newline <- rbind(newline, newline2)
  if (i==1) (genbl <- newline) else (genbl <- rbind(genbl, newline))
}
# new dataset: number of localities and genera for each region and global time bin
summary(genbl)
unique(genbl$Region)
summary(as.factor(genbl$GlobalTimeBins))
summary(as.factor(genbl$ContinentalBins))
# sort out zeros and NAs
genbl$GlobalTimeBins[which(genbl$GlobalTimeBins=="999")] <- NA
genbl$ContinentalBins[which(genbl$ContinentalBins=="999")] <- NA
for (i in 4:length(genbl)) (genbl[which(is.na(genbl[,i])),i] <- 0)
genbl[which(genbl$NumLocs==0),5] <- NA
# append the missing time and region information
genbl$GlobalBinStart <- ds$GlobalBinStart[match(genbl$GlobalTimeBins, ds$GlobalTimeBins)]
genbl$GlobalBinEnd <- ds$GlobalBinEnd[match(genbl$GlobalTimeBins, ds$GlobalTimeBins)]
genbl$GlobalBinCenter <- genbl$GlobalBinStart - (genbl$GlobalBinStart - genbl$GlobalBinEnd)/2
genbl$ContinBinStart <- ds$ContinBinStart[match(genbl$ContinentalBins, ds$ContinentalBins)]
genbl$ContinBinEnd <- ds$ContinBinEnd[match(genbl$ContinentalBins, ds$ContinentalBins)]
genbl$ContinBinCenter <- genbl$ContinBinStart - (genbl$ContinBinStart - genbl$ContinBinEnd)/2
genbl$Continent <- ds$Continent[match(genbl$Region, ds$Region)]
genbl[1:15,]
# get rid of the wrong continental bins for each continent
all(genbl$NumLocs[which(is.element(genbl$ContinentalBins, unicotime[1:16]) & genbl$Continent=="North America")]==0)
genbl <- genbl[-which(is.element(genbl$ContinentalBins, unicotime[1:16]) & genbl$Continent=="North America"),]
all(genbl$NumLocs[which(is.element(genbl$ContinentalBins, unicotime[17:length(unicotime)]) & genbl$Continent!="North America")]==0)
genbl <- genbl[-which(is.element(genbl$ContinentalBins, unicotime[17:length(unicotime)]) & genbl$Continent!="North America"),]
genbl[1:25,]
write.table(genbl, "Data_Regions.txt", row.names=FALSE, col.names=TRUE, sep="\t")


# genus dataset for the continents, with global and continental time bins

unireg <- unique(ds$Continent) # the continents
for (i in seq(along=unireg)) { # loop through the continents
  # global time bins
  dsi <- ds[which(ds$Continent==unireg[i]),] # dataset for region i
  numgen <- tapply(dsi$Genus, dsi$GlobalTimeBins, luni) # total number of genera within each global time bin in continent i
  numloc <- tapply(dsi$latlong, dsi$GlobalTimeBins, luni) # total number of localities within each global time bin in continent i
  # combine the data together to a new dataset
  newline <- data.frame(Continent=rep(unireg[i], length(unigltime)), GlobalTimeBins=unigltime, ContinentalBins=rep(999, length(unigltime)), NumLocs=numloc[match(unigltime, names(numloc))], NumGenera=numgen[match(unigltime, names(numgen))])
  newline$GlobalTimeBins <- as.character(newline$GlobalTimeBins)
  # continental time bins
  dsi <- ds[which(ds$Continent==unireg[i] & !is.na(ds$ContinentalBins)),] # dataset for continent i
  numgen <- tapply(dsi$Genus, dsi$ContinentalBins, luni) # total number of genera within each continental time bin in continent i
  numloc <- tapply(dsi$latlong, dsi$ContinentalBins, luni) # total number of localities within each continental time bin in continent i
  # combine the data together to a new dataset
  newline2 <- data.frame(Continent=rep(unireg[i], length(unicotime)), GlobalTimeBins=rep(999, length(unicotime)), ContinentalBins=unicotime, NumLocs=numloc[match(unicotime, names(numloc))], NumGenera=numgen[match(unicotime, names(numgen))])
  newline2$ContinentalBins <- as.character(newline2$ContinentalBins)
  newline <- rbind(newline, newline2)
  if (i==1) (genbl <- newline) else (genbl <- rbind(genbl, newline))
}
# new dataset: number of localities and genera for each region and global time bin
summary(genbl)
unique(genbl$Continent)
summary(as.factor(genbl$GlobalTimeBins))
summary(as.factor(genbl$ContinentalBins))
# sort out zeros and NAs
genbl$GlobalTimeBins[which(genbl$GlobalTimeBins=="999")] <- NA
genbl$ContinentalBins[which(genbl$ContinentalBins=="999")] <- NA
for (i in 4:length(genbl)) (genbl[which(is.na(genbl[,i])),i] <- 0)
genbl[which(genbl$NumLocs==0),5] <- NA
# append the missing time and region information
genbl$GlobalBinStart <- ds$GlobalBinStart[match(genbl$GlobalTimeBins, ds$GlobalTimeBins)]
genbl$GlobalBinEnd <- ds$GlobalBinEnd[match(genbl$GlobalTimeBins, ds$GlobalTimeBins)]
genbl$GlobalBinCenter <- genbl$GlobalBinStart - (genbl$GlobalBinStart - genbl$GlobalBinEnd)/2
genbl$ContinBinStart <- ds$ContinBinStart[match(genbl$ContinentalBins, ds$ContinentalBins)]
genbl$ContinBinEnd <- ds$ContinBinEnd[match(genbl$ContinentalBins, ds$ContinentalBins)]
genbl$ContinBinCenter <- genbl$ContinBinStart - (genbl$ContinBinStart - genbl$ContinBinEnd)/2
genbl[1:15,]
# get rid of the wrong continental bins for each continent
all(genbl$NumLocs[which(is.element(genbl$ContinentalBins, unicotime[1:16]) & genbl$Continent=="North America")]==0)
genbl <- genbl[-which(is.element(genbl$ContinentalBins, unicotime[1:16]) & genbl$Continent=="North America"),]
all(genbl$NumLocs[which(is.element(genbl$ContinentalBins, unicotime[17:length(unicotime)]) & genbl$Continent!="North America")]==0)
genbl <- genbl[-which(is.element(genbl$ContinentalBins, unicotime[17:length(unicotime)]) & genbl$Continent!="North America"),]
genbl[1:25,]
write.table(genbl, "Data_Continents.txt", row.names=FALSE, col.names=TRUE, sep="\t")


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

# Make community matrices with abundances for each genus (columns).
# -- in each region for each global time slice (rows)
# -- in each continent for each global time slice (rows)
# -- in each region for each continental time slice (rows)
# -- in each continent for each continental time slice (rows)
# At the same time, make presence-absence matrices for each genus
# in each locality.

# Then use the estimateR() quantitative model to estimate total richness
# from the abundance matrices.
# Also use the locality presence-absence matrices to estimate the species
# pool size.

# Unite all data to one regional and one continental dataset.

# set your working directory!
setwd("D://Data//Gammadiv//FilesForUpload")
# NOTE: the following code will overwrite some datasets included in the upload.
# To avoid this, change file names.

ds <- read.delim("Data_GenusOccurrence.txt", as.is=TRUE)
head(ds)
any(duplicated(ds))
# add some locality identifiers to the data
ds$latlong <- paste(ds$LAT, ds$LONG, sep="_") # unique locations
ds$LocCoTime <- paste(ds$ContinentalBins, ds$latlong, sep="_") # unique localities in land mammal ages: location-by-continental-time combination
ds$LocGlTime <- paste(ds$GlobalTimeBins, ds$latlong, sep="_") # unique localities in stratigraphic stages: location-by-global-time combination
# check the location numbers are all right
length(unique(ds$latlong))
luni <- function(x) (length(unique(x))) # function to count unique entries in a vector
tapply(ds$latlong, ds$Continent, luni) # how many localities are unique for each continent?

library(vegan)


# - STEP 1: make abundance and presence matrices

# --- regions, global time bins (stratigraphic stages)

# first, need a presence-absence matrix for each genus (columns)
# in each locality for each global time bin
is.present <- function(x) { # a function to check whether a genus is present in the dataset
  if (length(x)>0) (return(1)) else (return(length(x)))
}
ds$Genus <- as.factor(ds$Genus)
str(ds$Genus) # 690 genera
uniloc <- sort(unique(ds$LocGlTime)) # 1704 localities in the stratigraphic stages
for (i in seq(along=uniloc)) { # loop through the localities
  print(i)
  dsi <- ds[which(ds$LocGlTime==uniloc[i]),] # occurrences for the locality i
  mymati <- t(as.matrix(tapply(dsi$Genus, dsi$Genus, is.present))) # make the matrix for locality i
  rownames(mymati) <- uniloc[i]
  if (i==1) (mymat <- mymati) else (mymat <- rbind(mymat, mymati))
}
str(mymat)
occmatgl <- mymat

# then, turn this into an abundance matrix for each region
# where abundance is the occupancy of genera (number of localities in the region)
ds$RegionGlTime <- paste(ds$Region, ds$GlobalTimeBins, sep="_") 
unireg <- sort(unique(ds$RegionGlTime)) # 46 region-global time bin-combinations
for (i in seq(along=unireg)) {
  # the presence matrix (genus by locality) for just region i
  mati <- mymat[which(is.element(rownames(mymat), ds$LocGlTime[ds$RegionGlTime==unireg[i]])),]
  # get the abundances for just region i
  if (is.matrix(mati)) { # if there is more than one locality
    abmati <- t(as.matrix(colSums(mati, na.rm=TRUE)))
  } else { # if there is only one locality
    abmati <- t(as.matrix(as.numeric(!is.na(mati))))
    colnames(abmati) <- names(mati)
  }
  rownames(abmati) <- unireg[i]
  if (i==1) (abmat <- abmati) else (abmat <- rbind(abmat, abmati))
}
str(abmat)
abmat[1:10,1:20]
amreggl <- abmat


# --- continents, global time bins

# turn the presence matrix generated above into an abundance matrix for each continent
ds$ContGlTime <- paste(ds$Continent, ds$GlobalTimeBins, sep="_")
unireg <- sort(unique(ds$ContGlTime)) # 16 continent - global time bin - combinations
for (i in seq(along=unireg)) {
  # the presence matrix (genus by locality) for just region i
  mati <- mymat[which(is.element(rownames(mymat), ds$LocGlTime[ds$ContGlTime==unireg[i]])),]
  # get the abundances for just region i
  if (is.matrix(mati)) {
    abmati <- t(as.matrix(colSums(mati, na.rm=TRUE)))
  } else {
    abmati <- t(as.matrix(as.numeric(!is.na(mati))))
    colnames(abmati) <- names(mati)
  }
  rownames(abmati) <- unireg[i]
  if (i==1) (abmat <- abmati) else (abmat <- rbind(abmat, abmati))
}
str(abmat)
abmat[1:10,1:20]
amcontgl <- abmat


# --- regions, continental time bins (land mammal ages)

# apply the same code as above but for the smaller time bins.
# first, need a presence-absence matrix for each genus (columns)
# in each locality for each continental time bin
uniloc <- sort(unique(ds$LocCoTime)) # 1775 localities
for (i in seq(along=uniloc)) { # loop through the localities
  print(i)
  dsi <- ds[which(ds$LocCoTime==uniloc[i]),] # occurrences for the locality i
  mymati <- t(as.matrix(tapply(dsi$Genus, dsi$Genus, is.present))) # make the matrix for locality i
  rownames(mymati) <- uniloc[i]
  if (i==1) (mymat <- mymati) else (mymat <- rbind(mymat, mymati))
}
str(mymat)
head(mymat)
occmatco <- mymat

# then, turn this into an abundance matrix for each region
ds$RegionCoTime <- paste(ds$Region, ds$ContinentalBins, sep="_")
unireg <- sort(unique(ds$RegionCoTime)) # 89 region - continental time bin - combinations
unireg <- unireg[-grep("_NA", unireg)] # get rid of 3 - the North American data that could not be assigned a land mammal age (but was included above because it fell into a unique stratigraphic stage)
for (i in seq(along=unireg)) {
  # the presence matrix (genus by locality) for just region i
  mati <- mymat[which(is.element(rownames(mymat), ds$LocCoTime[ds$RegionCoTime==unireg[i]])),]
  # get the abundances for just region i
  if (is.matrix(mati)) {
    abmati <- t(as.matrix(colSums(mati, na.rm=TRUE)))
  } else {
    abmati <- t(as.matrix(as.numeric(!is.na(mati))))
    colnames(abmati) <- names(mati)
  }
  rownames(abmati) <- unireg[i]
  if (i==1) (abmat <- abmati) else (abmat <- rbind(abmat, abmati))
}
str(abmat)
abmat[1:10,1:20]
amregco <- abmat


# --- continents, continental time bins

# turn the presence matrix generated above into an abundance matrix for each continent
ds$ContCoTime <- paste(ds$Continent, ds$ContinentalBins, sep="_")
unireg <- sort(unique(ds$ContCoTime)) # 32 continent - continental time bin - combinations
unireg <- unireg[-grep("_NA", unireg)] # get rid of 1 - the North American data that does not fall into a unique time bin
for (i in seq(along=unireg)) {
  # the presence matrix (genus by locality) for just region i
  mati <- mymat[which(is.element(rownames(mymat), ds$LocCoTime[ds$ContCoTime==unireg[i]])),]
  # get the abundances for just region i
  if (is.matrix(mati)) {
    abmati <- t(as.matrix(colSums(mati, na.rm=TRUE)))
  } else {
    abmati <- t(as.matrix(as.numeric(!is.na(mati))))
    colnames(abmati) <- names(mati)
  }
  rownames(abmati) <- unireg[i]
  if (i==1) (abmat <- abmati) else (abmat <- rbind(abmat, abmati))
}
str(abmat)
abmat[1:10,1:20]
amcontco <- abmat


# save the occurrence and abundance matrices
save(occmatgl, occmatco, file="Output_GenusPresenceAbsenceMatrices.RData")
save(amreggl, amcontgl, amregco, amcontco, file="Output_GenusAbundanceMatrices.RData")

# to load these data back into R, use:
# load("Output_GenusPresenceAbsenceMatrices.RData")
# load("Output_GenusAbundanceMatrices.RData")


# - STEP 2: apply richness estimators
# get the species pool estimates and save them with the raw genus richness datasets 

# -- regional dataset

regds <- read.delim("Data_Regions.txt", as.is=TRUE)
# NOTE: if you ran  through the above script on compiling raw diversity data, Data_Regions as saved by you above will be the right dataset.
# If you start from here, Data_Regions from the uploaded filed should be cut to only the first 12 columns, otherwise you will duplicate columns that are already there.
# Add identifiers to match regional data to the matrices generated:
regds$RegionGlTime <- paste(regds$Region, regds$GlobalTimeBins, sep="_")
regds$RegionCoTime <- paste(regds$Region, regds$ContinentalBins, sep="_")
regds$RowNumber <- seq(along=regds[,1])
head(regds)
zeroize <- function (x) {  # a function to turn NAs into zeros in a matrix
  x[which(is.na(x))] <- 0
  return(x)
}

# --- regional dataset, global time bins

mymat2 <- apply(occmatgl, 2, zeroize) # turn NA into 0 in the presence-absence matrix
mymat2[1:10,1:10]
regpool <- ds$RegionGlTime[match(rownames(mymat2), ds$LocGlTime)] # classification of regions for the presence-absence matrix localities
head(regpool)
# from the occurrence matrix, estimate species pool size for the regions
poolreggl <- specpool(mymat2, pool=regpool)
# from the abundance matrix, estimate species pool size for the regions
estreggl <- t(estimateR(amreggl))
# merge the two datasets
reggl <- merge(estreggl, poolreggl[,-1], by=0, all=TRUE)
row.names(reggl) <- reggl$Row.names
reggl <- reggl[,-1]
# tidy up the the NAs and NaNs
if (any(is.na(reggl[,length(reggl)]))) (reggl[which(is.na(reggl[,length(reggl)])),length(reggl)] <- 0)
for (k in c(2:5,7,9,12)) {
  if (is.element(k, c(3:5))) {
    if (any(is.na(reggl[,k]))) (reggl[which(is.na(reggl[k])),k] <- NA)
  }
  if (any(reggl[which(!is.na(reggl[,k])),k]==0)) (reggl[which(!is.na(reggl[,k]) & reggl[,k]==0),k] <- NA)
}
head(reggl)
summary(reggl)
# merge with the regional dataset
regds2 <- merge(regds, reggl, by.x="RegionGlTime", by.y=0, all.x=TRUE, all.y=FALSE, sort=FALSE)
head(regds2)
all(regds2$NumLocs==regds2$n)

# --- regional dataset, continental time bins

mymat2 <- apply(occmatco, 2, zeroize) # turn NA into 0 in the presence-absence matrix
mymat2[1:10,1:10]
regpool <- ds$RegionCoTime[match(rownames(mymat2), ds$LocCoTime)] # classification of regions for the presence-absence matrix localities
head(regpool)
# from the occurrence matrix, estimate species pool size for the regions
poolregco <- specpool(mymat2, pool=regpool)
poolregco <- poolregco[-grep("_NA", rownames(poolregco)),] # get rid of the weird North American data that does not fall into a unique time bin
# from the abundance matrix, estimate species pool size for the regions
estregco <- t(estimateR(amregco))
# merge the two datasets
regco <- merge(estregco, poolregco[,-1], by=0, all=TRUE)
row.names(regco) <- regco$Row.names
regco <- regco[,-1]
# tidy up the the NAs and NaNs
if (any(is.na(regco[,length(regco)]))) (regco[which(is.na(regco[,length(regco)])),length(regco)] <- 0)
for (k in c(2:5,7,9,12)) {
  if (is.element(k, c(3:5))) {
    if (any(is.na(regco[,k]))) (regco[which(is.na(regco[k])),k] <- NA)
  }
  if (any(regco[which(!is.na(regco[,k])),k]==0)) (regco[which(!is.na(regco[,k]) & regco[,k]==0),k] <- NA)
}
summary(regco)
# merge with the regional dataset
for (i in seq(along=regco)) { # loop through the columns because this is a partial match process
  regds2[which(is.element(regds2$RegionCoTime, rownames(regco))),which(names(regds2)==names(regco)[i])] <- regco[match(regds2$RegionCoTime[which(is.element(regds2$RegionCoTime, rownames(regco)))], rownames(regco)),i]
}

# check this worked...
head(regds2)
names(regds2)
regds2[1,]
reggl["Western Europe_Aquitanian",]
regds2[which(regds2$RegionCoTime=="Western Europe_MN16"),]
regco["Western Europe_MN16",]
# sort out NAs for regions without localities
regds2$NumLocs[which(is.na(regds2$NumLocs))] <- 0    
regds2$n[which(is.na(regds2$n))] <- 0    
all(regds2$NumLocs==regds2$n)
all(regds2$NumGenera[!is.na(regds2$NumGenera)]==regds2$S.obs[!is.na(regds2$NumGenera)])
# sort out the dataset and save
regds2 <- regds2[order(regds2$RowNumber),]
names(regds2)
regds3 <- regds2[,c(2:13,17:27)]
head(regds3)
# NOTE: next line overwrites the dataset with the new one!
write.table(regds3, "Data_Regions.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# -- continental dataset

regds <- read.delim("Data_Continents.txt", as.is=TRUE)
# NOTE: if you ran  through the above script on compiling raw diversity data, Data_Continents as saved by you above will be the right dataset.
# If you start from here, Data_Continents from the upload should be cut to only the first 12 columns, otherwise you will duplicate columns that are already there.
# Add identifiers to match regional data to the matrices generated:
regds$RegionGlTime <- paste(regds$Continent, regds$GlobalTimeBins, sep="_")
regds$RegionCoTime <- paste(regds$Continent, regds$ContinentalBins, sep="_")
regds$RowNumber <- seq(along=regds[,1])
head(regds)

# --- continental dataset, global time bins

mymat2 <- apply(occmatgl, 2, zeroize) # turn NA into 0 in the presence-absence matrix
mymat2[1:10,1:10]
regpool <- ds$ContGlTime[match(rownames(mymat2), ds$LocGlTime)] # classification of regions for the presence-absence matrix localities
head(regpool)
# from the occurrence matrix, estimate species pool size for the continents
poolreggl <- specpool(mymat2, pool=regpool)
# from the abundance matrix, estimate species pool size for the continents
estreggl <- t(estimateR(amcontgl))
# merge the two datasets
reggl <- merge(estreggl, poolreggl[,-1], by=0, all=TRUE)
row.names(reggl) <- reggl$Row.names
reggl <- reggl[,-1]
# tidy up the the NAs and NaNs
if (any(is.na(reggl[,length(reggl)]))) (reggl[which(is.na(reggl[,length(reggl)])),length(reggl)] <- 0)
for (k in c(2:5,7,9,12)) {
  if (is.element(k, c(3:5))) {
    if (any(is.na(reggl[,k]))) (reggl[which(is.na(reggl[k])),k] <- NA)
  }
  if (any(reggl[which(!is.na(reggl[,k])),k]==0)) (reggl[which(!is.na(reggl[,k]) & reggl[,k]==0),k] <- NA)
}
summary(reggl)
# merge with the regional dataset
regds2 <- merge(regds, reggl, by.x="RegionGlTime", by.y=0, all.x=TRUE, all.y=FALSE, sort=FALSE)
head(regds2)

# --- continental dataset, continental time bins

mymat2 <- apply(occmatco, 2, zeroize) # turn NA into 0 in the presence-absence matrix
mymat2[1:10,1:10]
regpool <- ds$ContCoTime[match(rownames(mymat2), ds$LocCoTime)] # classification of regions for the presence-absence matrix localities
head(regpool)
# from the occurrence matrix, estimate species pool size for the continents
poolregco <- specpool(mymat2, pool=regpool)
poolregco <- poolregco[-grep("_NA", rownames(poolregco)),] # get rid of the weird North American data that does not fall into a unique time bin
# from the abundance matrix, estimate species pool size for the continents
estregco <- t(estimateR(amcontco))
# merge the two datasets
regco <- merge(estregco, poolregco[,-1], by=0, all=TRUE)
row.names(regco) <- regco$Row.names
regco <- regco[,-1]
# tidy up the the NAs and NaNs
if (any(is.na(regco[,length(regco)]))) (regco[which(is.na(regco[,length(regco)])),length(regco)] <- 0)
for (k in c(2:5,7,9,12)) {
  if (is.element(k, c(3:5))) {
    if (any(is.na(regco[,k]))) (regco[which(is.na(regco[k])),k] <- NA)
  }
  if (any(regco[which(!is.na(regco[,k])),k]==0)) (regco[which(!is.na(regco[,k]) & regco[,k]==0),k] <- NA)
}
summary(regco)
# merge with the continental dataset
for (i in seq(along=regco)) { # loop through the columns because this is a partial match process
  regds2[which(is.element(regds2$RegionCoTime, rownames(regco))),which(names(regds2)==names(regco)[i])] <- regco[match(regds2$RegionCoTime[which(is.element(regds2$RegionCoTime, rownames(regco)))], rownames(regco)),i]
}

# check this worked...
head(regds2)
summary(regds2)
regds2[1,]
reggl["Europe_Aquitanian",]
regds2[which(regds2$RegionCoTime=="Europe_MN16"),]
regco["Europe_MN16",]
# no continents without localities
regds2$NumLocs[which(is.na(regds2$NumLocs))]
regds2$n[which(is.na(regds2$n))]
all(regds2$NumLocs==regds2$n)
all(regds2$NumGenera[!is.na(regds2$NumGenera)]==regds2$S.obs[!is.na(regds2$NumGenera)])
# sort out the dataset and save
regds2 <- regds2[order(regds2$RowNumber),]
names(regds2)
regds3 <- regds2[,c(2:12,16:26)]
head(regds3)
# NOTE: next line overwrites the dataset with the new one!
write.table(regds3, "Data_Continents.txt", sep="\t", row.names=FALSE, col.names=TRUE)


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


# get the size of the area covered by mammalian localities
# for each continent and region in each global time bin (stratigraphic stage)

# At the same time, append the grid cell IDs for each
# locality to the Data_GenusOccurrence file.

# set your working directory!
setwd("D://Data//CatherineJanJussiProject//Gammadiv//FilesForUpload")
# NOTE: the following code will overwrite some datasets included in the upload.
# To avoid this, change file names.

ds <- read.delim("Data_GenusOccurrence.txt", as.is=TRUE)
head(ds)
any(duplicated(ds))
# add some locality identifiers to the data
ds$latlong <- paste(ds$LAT, ds$LONG, sep="_") # unique locations
# check the location numbers are all right
length(unique(ds$latlong))
luni <- function(x) (length(unique(x))) # function to count unique entries in a vector
tapply(ds$latlong, ds$Continent, luni) # how many localities are unique for each continent?
library(sp)
library(maps)
library(rgeos)
library(maptools)
myregs <- readShapePoly("Data_RegionsGridded.shp")
# AGAIN, NOTE: this dataset already contains the columns generated by the following code.
# If you want to rerun this, reduce the dataset to the first 9 columns:
# myregs@data <- myregs@data[,1:9]
names(myregs@data)
summary(as.factor(myregs@data$Region))
plot(myregs, border="grey")
map(add=TRUE)


# genus dataset for the regions, with global and continental time bins
# for each region, make a hull shape in each time bin from which area can then be calculated
# at the same time, generate a gridded dataset for the localities

unireg <- unique(ds$Region) # the regions in the dataset
unigltime <- unique(ds$GlobalTimeBins) # the stratigraphic stages
ds$GlobalTimeBins <- as.factor(ds$GlobalTimeBins)
fidcount <- 1 # for the FIDs, we need a counter that goes up after making each shape
for (i in seq(along=unireg)) { # loop through the regions
  print(paste(i))
  # convex hull for all localities in region i
  dsi <- ds[which(ds$Region==unireg[i]),c(4:5,10)] # dataset for region i: lat, long, global bin
  locsi <- SpatialPoints(dsi[!duplicated(dsi[,1:2]),c(2,1)]) # locations in regions i through all time bins
  reghull <- gConvexHull(locsi) # get the polygon containing all sampling locations
  mypols <- SpatialPolygonsDataFrame(reghull, data=data.frame(Region=unireg[i], GlobalTimeBins="All", COntinentalBins="All", stringsAsFactors=FALSE))
  if (i>1) {
    mypols <- spChFIDs(mypols, as.character(fidcount)) # fiddle the FIDs
    pols <- rbind(pols, mypols) # make the polygon dataset to hold all polygons
  } else (pols <- mypols)
  fidcount <- fidcount + 1 # keep track of the FID to assign
  # hulls and grid datasets for each of the global time bins
  for (k in seq(along=unigltime)) {
    locsik <- dsi[which(dsi$GlobalTimeBins==unigltime[k]),c(2,1)] # the locations in global time slice k
    if (length(locsik[!duplicated(locsik),1])>0) { # spatial points cannot be made with empty datasets
      pointsik <- SpatialPoints(locsik[!duplicated(locsik),]) # localities in region i and time bin k
      if (length(locsik[!duplicated(locsik),1])>2) { # if there are more than 2 localities (can't make a polygon otherwise)
        # make the polygon around the mammalian localities and save it
        polik <- SpatialPolygonsDataFrame(gConvexHull(pointsik), data=data.frame(Region=unireg[i], GlobalTimeBins=unigltime[k], stringsAsFactors=FALSE))
        polik <- spChFIDs(polik, as.character(fidcount)) # fiddle the FIDs
        fidcount <- fidcount + 1 # keep track of the FID to assign
        pols <- rbind(pols, polik) # add the hull for bin k to the polygon dataset
      }
      # get the information with location is in which 1 degree grid cell
      try1 <- over(pointsik, myregs, returnList=FALSE)
      if (any(as.character(try1$Region)!=unireg[i])) { # if a point falls on the edge between regions...
        print(paste(locsik[which(as.character(try1$Region)!=unireg[i]),])) # ...print out which one it was so I can fix it later
      }
      # if all is fine, no printing occurs and the dataset gets generated
      griddsik <- cbind(locsik[!duplicated(locsik),], try1) # a dataset with locality lat and long and the information which grid cell this falls into
      griddsik$GlobalTimeBin <- unigltime[k] # add the time bin k
      if (all(i==1, k==1)) (gridds <- griddsik) else (gridds <- rbind(gridds, griddsik))
    }
  }
}

# check the grid-locality dataset
head(gridds)
table(gridds$Region, gridds$GlobalTimeBin) # check number of localities
# append the grid cell ID information to the occurrence dataset
ds$latlongbin <- paste(ds$latlong, ds$GlobalTimeBins, sep="_")
gridds$latlongbin <- paste(gridds$LAT, gridds$LONG, gridds$GlobalTimeBin, sep="_")
all(is.element(ds$latlongbin, gridds$latlongbin))
all(is.element(gridds$latlongbin, ds$latlongbin))
ds$WM_ID <- gridds$WM_ID[match(ds$latlongbin, gridds$latlongbin)]
ds$Centroid_X <- gridds$Centroid_X[match(ds$latlongbin, gridds$latlongbin)]
ds$Centroid_Y <- gridds$Centroid_Y[match(ds$latlongbin, gridds$latlongbin)]
ds$Region2 <- gridds$Region[match(ds$latlongbin, gridds$latlongbin)]
all(ds$Region==ds$Region2) # just checking it assigned the same regions
names(ds)
ds2 <- ds[,c(1:20,23:25)]
# NOTE: next line overwrites the dataset with the new one!
write.table(ds2, "Data_GenusOccurrence.txt", sep="\t", row.names=FALSE, col.names=TRUE)

# check and save the hulls
str(pols@data)
plot(subset(pols, pols@data$GlobalTimeBins=="All"))
map(add=TRUE, col="grey")
plot(subset(pols, pols@data$GlobalTimeBins!="All"), add=TRUE, border="red")
summary(as.factor(pols$GlobalTimeBins))
summary(as.factor(pols$COntinentalBins))
writePolyShape(pols, fn="RegionsHulls")

# I edited this file in ArcGIs:
# - Western and Eastern European polygons sometimes overlapped, this got fixed
# - calculated area size in ArcGIS with Mollweide projection
# the edited and projected data Data_RegionsHulls_edproj is included
# in the uploaded files.


# append the area information from the hulls to the regional and continental datasets

# regional dataset: load the data so far
ds <- read.delim("Data_Regions.txt", as.is=TRUE)
head(ds)
# generate empty columns for the area data
ds$Area <- 9999999
ds$Area2 <- 9999999
library(foreign)
# load the edited hulls with area column
ars <- read.dbf("RegionsHulls_edproj.dbf")
head(ars)
table(ars$Region, ars$GlobalTime) # check which regions don't have a hull
ars <- ars[which(ars$GlobalTime!="All"),] # the hull across all time bins is not needed
# modify the area of western North America
# because of the tectonic extension
ars$Area2 <- ars$Area
ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Aquitanian"] <- ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Aquitanian"] * 0.5
ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Burdigalian"] <- ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Burdigalian"] * 0.5
ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Langhian"] <- ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Langhian"] * 0.6
ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Serravalian"] <- ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Serravalian"] * 0.7
ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Tortonian"] <- ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Tortonian"] * 0.8
ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Messinian"] <- ars$Area2[ars$Region=="Western North America" & ars$GlobalTime=="Messinian"] * 0.9
unireg <- unique(ds$Region) # loop through the unique regions
for (i in seq(along=unireg)) {
  arsi1 <- ars[which(ars$Region==unireg[i]),] # area dataset for region i
  # put the right area values into the column in the region dataset
  ds$Area[which(ds$Region==unireg[i] & is.element(ds$GlobalTimeBins, arsi1$GlobalTime))] <- arsi1$Area[match(ds$GlobalTimeBins[which(ds$Region==unireg[i] & is.element(ds$GlobalTimeBins, arsi1$GlobalTime))], arsi1$GlobalTime)]
  ds$Area2[which(ds$Region==unireg[i] & is.element(ds$GlobalTimeBins, arsi1$GlobalTime))] <- arsi1$Area2[match(ds$GlobalTimeBins[which(ds$Region==unireg[i] & is.element(ds$GlobalTimeBins, arsi1$GlobalTime))], arsi1$GlobalTime)]
}
summary(ds$Area) # have to set the NAs
ds$Area[which(ds$Area==9999999)] <- NA
ds$Area2[which(ds$Area2==9999999)] <- NA
summary(ds$Area) # lots of NAs because the dataset contains the continental time bins too, for which I did not calculate area
summary(ds$Area2)
summary(ds$Area[!is.na(ds$GlobalTimeBins)]) # this should match the numbers expected from the table across "ars" dataset above
summary(ds$Area2[!is.na(ds$GlobalTimeBins)])
head(ds)
# NOTE: next line overwrites the dataset with the new one!
write.table(ds, "Data_Regions.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# continental data: generate this from the regional data
ars$Continent <- ds$Continent[match(ars$Region, ds$Region)] # get the continent column into the ars dataset from the ds dataset by matching the regions
head(ars)
summary(as.factor(ars$Continent))
# continental dataset: load the data so far
ds <- read.delim("Data_Continents.txt", as.is=TRUE)
head(ds)
# generate empty columns for the area data
ds$Area <- 9999999
ds$Area2 <- 9999999
unireg <- unique(ds$Continent) # loop through the unique continents
for (i in seq(along=unireg)) {
  arsi <- ars[which(ars$Continent==unireg[i]),] # area dataset for continent i
  myareas <- tapply(arsi$Area, arsi$GlobalTime, sum) # calculate the area as the sum across regions within the continent, for each global time bin
  # put the right area values into the column in the region dataset
  ds$Area[which(ds$Continent==unireg[i] & is.element(ds$GlobalTimeBins, names(myareas)))] <- myareas[match(ds$GlobalTimeBins[which(ds$Continent==unireg[i] & is.element(ds$GlobalTimeBins, names(myareas)))], names(myareas))]
  # same for area 2 (the adjusted values for tectonic extension of North America)
  myareas2 <- tapply(arsi$Area2, arsi$GlobalTime, sum)
  ds$Area2[which(ds$Continent==unireg[i] & is.element(ds$GlobalTimeBins, names(myareas2)))] <- myareas2[match(ds$GlobalTimeBins[which(ds$Continent==unireg[i] & is.element(ds$GlobalTimeBins, names(myareas2)))], names(myareas2))]
}
summary(ds$Area) # have to set the NAs
ds$Area[which(ds$Area==9999999)] <- NA
ds$Area2[which(ds$Area2==9999999)] <- NA
summary(ds$Area) # lots of NAs because the dataset contains the continental time bins too, for which I did not calculate area
summary(ds$Area2)
summary(ds$Area[!is.na(ds$GlobalTimeBins)]) # this should have no NAs
summary(ds$Area2[!is.na(ds$GlobalTimeBins)])
head(ds)
# NOTE: next line overwrites the dataset with the new one!
write.table(ds, "Data_Continents.txt", sep="\t", row.names=FALSE, col.names=TRUE)

