###
###  Read and process the OBIS bivalve data set
### 
###  Before this script, run the "worms_prepare.R" script.
###
###
## 1. Libraries, functions
library(robis)
library(finch)
library(dplyr)
mutate_cond <- function(.data, condition, ..., envir = parent.frame()) {
  condition <- eval(substitute(condition), .data, envir)
  .data[condition, ] <- .data[condition, ] %>% mutate(...)
  .data
}
#
# Load OBIS data:
# ob <- occurrence("Bivalvia")
# 
# read OBIS data when already saved. Otherwise download from https://obis.org/
#
setwd("E://Erlangen_2021/OBIS/data_and_code/data")
obiv <- readRDS("OBIS_bivalves_download_2020_10_26.rds")
#
# Family and species level mineralogy
bivmin <- read.csv("bivalve_families_min_and_infaunal.csv")
bivminspec <- read.csv("Bivalve_species_mineralogies.csv")
#
obisfam <- as.character(unique(obiv$family))
wormsfam <- bivmin$family
#
# correct families
# obiv$family[which(obiv$family == "Modiolidae")] <- "Mytilidae"
# obiv$family[which(obiv$family == "Montacutidae")] <- "Lasaeidae"
# obiv$family[which(obiv$family == "Kelliidae")] <- "Lasaeidae"
obiv$family[which(obiv$family == "Bathyspinulidae")] <- "Nuculanidae"
# obiv$family[which(obiv$family == "Crenellidae")] <- "Mytilidae"
# obiv$family[which(obiv$family == "Leptonidae")] <- "Lasaeidae"
#
# problematic families: Montacutidae (unaccepted, superfamily Galeommatoidea): 50 000
# problematic families: Kelliidae (unaccepted, superfamily Galeommatoidea):     1 000
# problematic families: Bathyspinulidae (unaccepted, = family Ledellinae):        300
# problematic families: Crennellidae (unaccepted, superfamily Mytiloidea):          3
# problematic families: Manzanellidae (extinct):                                   53
# problematic families: Mactromyidae (extinct):                                    15
# problematic families: Inoceramida (extinct):                                     11
# problematic families: Leptonidae (unaccepted, = family Lasaeidae):               10
# problematic families: Schizodidae (extinct):                                      7
# problematic families: Radiolitidae (extinct):                                     6
# problematic families: Plagioptychidae (extinct):                                  1
# problematic families: Plagioptychidae (unaccepted, = subfamily Modiolinae, superfamily Mytiloidea): 1
#
# single out extinct families
extinctfam <- c("Manzanellidae", "Mactromyidae", "Inoceramida", "Schizodidae", "Radiolitidae", 
                "Plagioptychidae", "Halobiidae", "Iotrigoniidae", "Monotidae")
#
# 1339903 occurrences
# - 87 extinct occurrences
# remove extinct
obiv <- obiv %>% filter(!(family %in%  extinctfam))
#
# remove unresolved families (- 65388 occurrences)
obiv <- obiv %>% filter(!(family == ""))
#
# 1274428 remain
#########################################################
#
# assign mineralogy based on family
#
whatmin <- function(zaragonite,zcalcite) {
  if(zaragonite == 1 & zcalcite == 0) out <- "ara only" 
  if(zaragonite == 1 & zcalcite == 1) out <- "mixed"
  if(zaragonite == 0 & zcalcite == 1) out <- "cal only"
  out
}
#
fa_min_obiv <- rep(NA,nrow(obiv))
current_min <- NA
obiv$family <- as.character(obiv$family)
for (i in 1:nrow(bivmin)) {
  current_fam <- bivmin[i,]
  current_min <- whatmin(current_fam$aragonite, current_fam$calcite)
  fa_min_obiv[which(obiv$family == bivmin$family[i])] <- current_min
}
table(fa_min_obiv)
#
min_obiv <- fa_min_obiv
#
obiv$genus <- as.character(obiv$genus)
obiv$species <- as.character(obiv$species)
# assign mineralogy based on species
for (i in 1:(nrow(bivminspec)-1)) {
  current_spec <- bivminspec[i,]
  current_min <- whatmin(current_spec$aragonite, current_spec$calcite)
  min_obiv[which(obiv$species == paste(as.character(current_spec$genus), as.character(current_spec$species)))] <- current_min
}
#
table(min_obiv)
#
min_obiv[which(obiv$family == "Mytilidae" & 
                 !(obiv$species %in% paste(as.character(bivminspec$genus), as.character(bivminspec$species))))] <- "unknown"
#
obiv <- obiv %>% mutate(mineral = min_obiv)
#
obiv$mineral[which(obiv$genus=="Tridacna")] <- "ara only"
#
### infaunal ones 
#
fa_infaunal_obiv <- rep(NA,nrow(obiv))
current_infaunal <- NA
obiv$family <- as.character(obiv$family)
for (i in 1:nrow(bivmin)) {
  current_fam <- bivmin[i,]
  current_infaunal <-  as.numeric(as.character(current_fam$infaunal))
  fa_infaunal_obiv[which(obiv$family == current_fam$family)] <- current_infaunal
}
table(fa_infaunal_obiv)
#
obiv <- obiv %>% mutate(infaunal = fa_infaunal_obiv)
# Correct the genus Tridacna to be epifaunal
obiv$infaunal[which(obiv$genus == "Tridacna")] <- 0
#
table(is.na(obiv$species))
#
###########################################################################################
#
# Species level data set
#
obiv <- obiv %>% filter(!(species == ""))
# removes 136582 occurrences, 1137846 remain
#
#### OBIS bivalve species environmental range
#
allspecies <- unique(obiv$species)
sort(allspecies)
#
# check if accepted (needs to have the WORMS data loaded!)
worms_accepted_species <- (biv %>% filter(taxonRank == "species"))$acname 
worms_unaccepted_species <- (biv_all %>% filter(taxonRank == "species" & scientificName != acname)) 
#
# View((biv %>% filter(taxonRank == "species")))
# table(allspecies %in% worms_accepted_species)
wrongspecies <- allspecies[which(!(allspecies %in% worms_accepted_species))]
# sort(wrongspecies)
# table(wrongspecies %in% worms_unaccepted_species$scientificName)

# length(unlist(c(sapply(wrongspecies, function(x) worms_unaccepted_species$acname[which(worms_unaccepted_species$scientificName==x)]))))

### Correct species names
#
# Omit missing species occurrences
#obiv <- subset(obiv,!(is.na(species)))
#nrow(obiv)
#
# Correct the ones that are corrected in Worms

for(i in 1:length(wrongspecies)) {
  if(wrongspecies[i] %in% worms_unaccepted_species$scientificName) {
    obiv$species[which(obiv$species==wrongspecies[i])] <- rep(worms_unaccepted_species$acname[which(worms_unaccepted_species$scientificName==wrongspecies[i])], 
                                                            length(obiv$species[which(obiv$species==wrongspecies[i])]))
  }
  print(i)
}

# unique(obiv$species)[which(!(unique(obiv$species) %in% worms_accepted_species))] %in% allspecies 

still_missing <- wrongspecies[which(!(wrongspecies %in% worms_unaccepted_species$scientificName))]


obiv$species[which(obiv$species == "Paracratis minuta")] <- "Limopsis minuta"
obiv$species[which(obiv$species == "Paracratis natalis")] <- "Limopsis natalis"
obiv$species[which(obiv$species == "Oblimopa compressa")] <- "Limopsis compressa"
obiv$species[which(obiv$species == "Pecten benedictus")] <- "Spathochlamys benedicti"
obiv$species[which(obiv$species == "Aspalima affinis")] <- "Limopsis affinis"
obiv$species[which(obiv$species == "Parvamussium jeffreysi")] <- "Parvamussium jeffreysii"
obiv$species[which(obiv$species == "Aspalima cristata")] <- "Limopsis cristata"
obiv$species[which(obiv$species == "Aspalima mabilliana")] <- "Limopsis mabilliana"
obiv$species[which(obiv$species == "Chlamys (Zygochlamys) lischkei")] <- "Zygochlamys patagonica"
obiv$species[which(obiv$species == "Volachlamys tranquebarica")] <- "Volachlamys tranquebarica"  #### is accepted
obiv$species[which(obiv$species == "Pecten (Pseudamusium) sigsbeei")] <- "Pectinella sigsbei"
obiv$species[which(obiv$species == "Isognomon australicum")] <- "Isognomon australica"
obiv$species[which(obiv$species == "Parvamussium holmesi")] <- "Parvamussium holmesii"
obiv$species[which(obiv$species == "Ctenoides sanctipauli")] <- "Ctenoides sanctipauli"  #### is accepted
obiv$species[which(obiv$species == "Cardites antiquata")] <- "Cardites antiquatus"
obiv$species[which(obiv$species == "Anadara chemnitzi")] <- "Anadara chemnitzii"
obiv$species[which(obiv$species == "Limea (Isolimea) limopsis")] <- "Limea limopsis"
obiv$species[which(obiv$species == "Pycnodonte (Pycnodonte) taniguchii")] <- "Pycnodonte taniguchii"
obiv$species[which(obiv$species == "Aspalima galatheae")] <- "Limopsis galatheae"
obiv$species[which(obiv$species == "Hyotissa fisheri")] <- "Hyotissa fisheri" #### Taxon inquirendum
obiv$species[which(obiv$species == "Pecten speciosa")] <- "Gloripallium speciosum"
obiv$species[which(obiv$species == "Cardita floridana")] <- "Cardites floridanus"
obiv$species[which(obiv$species == "Hyalopecten parvulinus")] <- "Hyalopecten parvulinus" #### Taxon inquirendum
obiv$species[which(obiv$species == "Aspalima lanceolata")] <- "Limopsis lanceolata"
obiv$species[which(obiv$species == "Limea argentineae")] <- "Limea argentinae"
obiv$species[which(obiv$species == "Aspalima intermedia")] <- "Limopsis intermedia"
obiv$species[which(obiv$species == "Aspalima spicata")] <- "Limopsis spicata"
obiv$species[which(obiv$species == "Lima colorata")] <- "Lima zealandica" #### colorata is extinct according to worms, might be "colorata zealandica"
obiv$species[which(obiv$species == "Aspalima lata")] <- "Limopsis lata"
obiv$species[which(obiv$species == "Barbatia novaezelandiae")] <- "Barbatia novaezealandiae"
obiv$species[which(obiv$species == "Oblimopa woodwardi")] <- "Limopsis woodwardi"
obiv$species[which(obiv$species == "Oblimopa forteradiata")] <- "Limopsis forteradiata"
obiv$species[which(obiv$species == "Oblimopa tenuiradiata")] <- "Limopsis tenuiradiata"
obiv$species[which(obiv$species == "Barbatia coelata")] <- "Barbatia coelata" #### Taxon inquirendum
obiv$species[which(obiv$species == "Pinctada imbricata fucata")] <- "Pinctada fucata"
obiv$species[which(obiv$species == "Aspalima longipilosa")] <- "Limopsis longipilosa"
obiv$species[which(obiv$species == "Paracratis minuta")] <- "Limopsis minuta"
obiv$species[which(obiv$species == "Lampsilis ochracea")] <- "Atlanticoncha ochracea"

# sapply(still_missing, function(x) length(obiv$species[which(obiv$species == x)]))

unaccepted <- c("Poromya carinata"
                
)


correct_species <- 
  c("Gigantidas platifrons",
    "Periglypta exclathrata",
    "Papillicardium minimum",
    "Thyasira gouldii",
    "Phlyctiderma semiasperum",
    "Turneroconcha magnifica",
    "Gigantidas childressi",
    "Sarepta natalensis",
    "Tenuicorbula lyoni",
    "Archivesica kawamurai",
    "Mytella strigata",
    "Ostrea permollis",
    "Bathyspinula vityazi",
    "Nanostrea pinnicola",
    "Eucallista purpurata",
    "Rugalucina vietnamica",
    "Ostrea algoensis",
    "Notocina falklandica",
    "Phlyctiderma verrilli",
    "Solamen spectabile",
    "Serripecten tiorioriensis",
    "Cetoconcha spinosula",
    "Rugalucina angela",
    "Hawaiarca sulcata",
    "Papyridea clenchi",
    "Corbicula striatella",
    "Spissatella maxwelli",
    "Corbicula japonica",
    "Dallocardia delicatula",
    "Crassicardia crassidens",
    "Ctenoides sanctipauli",
    "Abditoconus heterosiphon",
    "Unio pictorum",
    "Joannisiella tsuchii",
    "Aligena ovalis",
    "Archivesica myriamae",
    "Solamen striatissimum",
    "Psammacoma inexpectata",
    "Archivesica nanshaensis",
    "Corbicula siamensis",
    "Papyridea curvicostata",
    "Varicorbula rotalis",
    "Lima waipipiensis",
    "Spondylus erectospinus",
    "Abditoconus anselli",
    "Gigantidas mauritanicus",
    "Myonera atlasiana",
    "Nivanteredo coronata",
    "Lithoredo abatanica",
    "Miltha neozelanica"
    )

extinct_species <- c("Amusium papyraceum",
                     "Yoldiella micrometrica",
                     "Panopea abrupta",
                     "Pecten benedictus",
                     "Serripecten yahliensis",
                     "Glycymerita manaiaensis",
                     "Mesopeplum costatostriatum",
                     "Neilonella pusio",
                     "Pteromyrtea dispar",
                     "Mancikellia divae",
                     "Ledella inflatella",
                     "Barytellina tertia",
                     "Mactra patagonica",
                     "Tucetona aucklandica",
                     "Finlayella imbellica",
                     "Pycnodonte tarda",
                     "Gari aucklandica",
                     "Lima colorata",
                     "Yoldiella ardgowanica",
                     "Lutraria trapezoidalis",
                     "Neilo wairoana",
                     "Semeloidea miocenica",
                     "Limatula trulla",
                     "Eucrassatella ampla",
                     "Panopea worthingtoni",
                     "Thyasira peruviana",
                     "Lentipecten hochstetteri",
                     "Pteromyrtea motutaraensis",
                     "Glycymerita huttoni",
                     "Thyasira marwicki",
                     "Barbatia gibba",
                     "Anadara lienosa",
                     "Spissatella trailli",
                     "Spisula crassitesta",
                     "Perrierina ovata",
                     "Spissatella clifdenensis",
                     "Crassostrea ingens",
                     "Thyasira motutaraensis",
                     "Spissatella poroleda",
                     "Serratina ferrari",
                     "Eucrassatella scopalveus",
                     "Glycymerita hurupiensis",
                     "Glycymeris shrimptoni",
                     "Mesopeplum kaiparaense",
                     "Lucina columbella",
                     "Pteria oneroaensis",
                     "Semeloidea donaciformis",
                     "Spissatella media",
                     "Glycymeris mahiana",
                     "Pycnodonte mackayi",
                     "Draculamya pumila",
                     "Finlayella sinuaris",
                     "Pinna distans",
                     "Flemingostrea wollastoni",
                     "Perna tetleyi",
                     "Chokekenia nicolasi",
                     "Glycymerita robusta",
                     "Myadora delta",
                     "Tucetona wairarapaensis",
                     "Megaxinus transversus",
                     "Athlopecten athleta",
                     "Crenostrea cannoni",
                     "Glycymeris waipipiensis",
                     "Lucinoma taylori",
                     "Parapholas thomsoni",
                     "Pteromyrtea laminata",
                     "Zemysia ampla",
                     "Thyasira minoensis",
                     "Myrtea microlirata",
                     "Mactra acuminella",
                     "Mesopeplum burnetti",
                     "Borniola kaawaensis",
                     "Gari oamarutica",
                     "Pseudoportlandia solenelloides",
                     "Barytellina crassidens",
                     "Glycymerita thomsoni",
                     "Solecurtus evolutus",
                     "Myllita finlayi",
                     "Parvamussium vafer",
                     "Arthritica finlayi",
                     "Glycymeris callaghani",
                     "Solecurtus bensoni",
                     "Crenostrea wuellerstorfi",
                     "Scalpomactra continua",
                     "Spisula voyi",
                     "Phialopecten triphooki",
                     "Spissatella acculta",
                     "Triplicitella australis",
                     "Saccella semiteres",
                     "Zygochlamys seymouri",
                     "Spondylus aucklandicus",
                     "Janupecten uttleyi",
                     "Zygochlamys geminata",
                     "Janupecten polemicus",
                     "Lutraria grandis",
                     "Mesopeplum crawfordi",
                     "Saccella probellula",
                     "Tucetona traversi",
                     "Gonimyrtea bucculenta",
                     "Thyasira planata",
                     "Notocallista watti",
                     "Bankia turneri",
                     "Mesopeplum waikohuense",
                     "Ctenoides naufragus",
                     "Saccella motutaraensis",
                     "Tellinota arta",
                     "Paphies otaiana",
                     "Zygochlamys jorgensis",
                     "Phialopecten marwicki",
                     "Glycymerita kaawaensis",
                     "Kurinuia areolata",
                     "Sectipecten allani",
                     "Tellinota tumens",
                     "Eucrassatella subobesa",
                     "Mesopeplum syagrus"
                     
                     )

taxon_inquirendum <- 
  c("Lyonsia adriatica",
    "Mysella subquadrata",
    "Nuculana fossa",
    "Magallana cuttackensis",
    "Pandora cistula",
    "Tellina barhampurensis",
    "Nuculana radiata",
    "Montacuta acuminata",
    "Pharella ovalis",
    "Magallana rivularis",
    "Gastrochaena indistincta",
    "Magallana gryphoides",
    "Semele isosceles",
    "Barbatia coelata",
    "Amussium solitarium",
    "Diplodonta annandalei",
    "Hyalopecten parvulinus",
    "Venericardia columnaria",
    "Hyotissa fisheri",
    "Astarte striata",
    "Donax nuxfagus",
    "Cumingia hinduorum",
    "Thyasira dilatata"
    
    )

nomen_dubium <- c("Chama japonica",
                  "Nucula turgida",
                  "Cardiomya knudseni",
                  "Psammotreta sanguinolenta",
                  "Tellina vernalis",
                  "Mysella minuscula",
                  "Chama fragum",
                  "Genaxinus bongraini",
                  "Tellina subangulata",
                  "Anodontia edentula",
                  "Erycina firmata",
                  "Modiola opifex",
                  "Carditopsis malvinae",
                  "Teredo heaphyi",
                  "Scacchia ovata",
                  "Psammotreta angulata",
                  "Cymatioa pulchra",
                  "Thyasira subtrigona",
                  "Lyonsia argentea",
                  "Loripes pertenuis",
                  "Diplodonta insulsa",
                  "Scrobicularia rubiginosa",
                  "Lucinoma antarctica",
                  "Solen acinaces"
                  )

not_found <- c("Limea tanzaniensis",
               "Limaria imitans",
               "Limatula surugensis",
               "Limaria pseudoorientalis",
               "Limaria burryi"
               
               )



temporary_name <- c("Vesicomya venusta",
                    "Kelliella pilula",
                    "Veneriglossa vesica",
                    "Kelliella concentrica"  
                    )

## remove extinct:
obiv <- obiv %>% filter(!(species %in% extinct_species))

## remove unaccepted:
obiv <- obiv %>% filter(!(species %in% unaccepted))

## remove not found:
obiv <- obiv %>% filter(!(species %in% not_found))

## remove nomen dubium:
obiv <- obiv %>% filter(!(species %in% nomen_dubium))

## remove taxon inquirendum:
obiv <- obiv %>% filter(!(species %in% taxon_inquirendum))

# unknown minerlaogy / not valid species
obiv <- subset(obiv, !(species %in% c("Chama fragum", "Carditopsis malvinae", "Chama japonica", "Limea tanzaniensis", "")))
obiv <- obiv %>% filter(originalScientificName != "Pecten nanus")


#################################################
#
# remove imprecise localities 

decimalplaces <- function(x) {
  if (abs(x - round(x)) > .Machine$double.eps^0.5) {
    nchar(strsplit(sub('0+$', '', as.character(x)), ".", fixed = TRUE)[[1]][[2]])
  } else {
    return(0)
  }
}
library(tidyr)
obiv <- obiv %>% mutate(latdecimals = sapply(obiv$decimalLatitude,function(x) decimalplaces(x)))
obiv$decimalLongitude[55938] <- 0.001
obiv$decimalLongitude[226870] <- 0.001
obiv$decimalLongitude[which(obiv$decimalLongitude < 0.001 & obiv$decimalLongitude > -0.001 & !(obiv$decimalLongitude==0))] <-
  0.001
obiv <- obiv %>% mutate(londecimals = sapply(obiv$decimalLongitude,function(x) decimalplaces(x)))

# exclude imprecise: - 25880 occurrences
obiv <- obiv[which(obiv$londecimals>=1 & obiv$latdecimals>=1),]


### Cleaning done

# remove infaunal

obiv <- subset(obiv,infaunal == 0)

# 279254 remain

# remove unknown mineralogy
obiv <- subset(obiv, mineral != "unknown" )

# 239276 remain


### Download Environmental data from the Bio-Oracle

library(sdmpredictors)
library(scales)
library(dplyr)

# Download specific layers to the current directory 
bt <- load_layers(c("BO2_tempmean_bdmean"))#, "BO2_tempmin_bdmean", "BO2_tempmax_bdmean")) 
# layers_correlation(c("BO2_tempmean_bdmean", "BO2_tempmin_bdmean", "BO2_tempmax_bdmean")) 

bs <- load_layers(c("BO2_salinitymean_bdmean"))#, "BO2_salinitymin_bdmean", "BO2_salinitymax_bdmean")) 
#layers_correlation(c("BO2_salinitymean_bdmean", "BO2_salinitymin_bdmean", "BO2_salinitymax_bdmean")) 

btr <- load_layers(c("BO2_temprange_bdmean")) 

btmin <- load_layers(c("BO2_tempmin_bdmean"))

btmax <- load_layers(c("BO2_tempmax_bdmean"))

obiv <- obiv %>% mutate(cell = (extract(bs, cellnumber = T, df = T, y = cbind(decimalLongitude,decimalLatitude)))$cells)

bathy_min <- load_layers(c("BO_bathymin"))

bathy_max <- load_layers(c("BO_bathymax"))

bathy_mean <- load_layers(c("BO_bathymean"))

sstmean_lay <- load_layers(c("BO21_tempmean_ss")) 
sstmin_lay <- load_layers(c("BO21_tempmin_ss")) 
sstmax_lay <- load_layers(c("BO21_tempmax_ss")) 
sstrange_lay <- load_layers(c("BO21_temprange_ss")) 

btr_min <- load_layers(c("BO2_temprange_bdmin")) 

btmin_min <- load_layers(c("BO2_tempmin_bdmin"))

btmax_min <- load_layers(c("BO2_tempmax_bdmin"))

bt_min <- load_layers(c("BO2_tempmean_bdmin"))

ss_sal <- load_layers(c("BO2_salinitymean_ss"))


# Test if the salinity and temperature grid cells give the same numbers
cbind((obiv[1:100,] %>% mutate(cell = (extract(bs, cellnumber = T, df = T, y = cbind(decimalLongitude,decimalLatitude)))$cells))$cell,
      (obiv[1:100,] %>% mutate(cell = (extract(bt, cellnumber = T, df = T, y = cbind(decimalLongitude,decimalLatitude)))$cells))$cell)
# yes

obiv <-  obiv %>% mutate_if(is.factor, as.character) # convert factors to characters

# Assign temperature and salinity to occurrences
obiv <-  obiv %>% mutate(salinity = raster::extract(bs,cell))
obiv <-  obiv %>% mutate(temperature = raster::extract(bt,cell))
obiv <-  obiv %>% mutate(temperaturerange = raster::extract(btr,cell))
obiv <-  obiv %>% mutate(temperatureminimum = raster::extract(btmin,cell))

obiv <-  obiv %>% mutate(BO_bathymin = raster::extract(bathy_min,cell))
obiv <-  obiv %>% mutate(BO_bathymax = raster::extract(bathy_max,cell))
obiv <-  obiv %>% mutate(BO_bathymean = raster::extract(bathy_mean,cell))

obiv <-  obiv %>% mutate(temperaturemin = raster::extract(btmin,cell))
obiv <-  obiv %>% mutate(temperaturemax = raster::extract(btmax,cell))

obiv <-  obiv %>% mutate(sstmin = raster::extract(sstmin_lay,cell))
obiv <-  obiv %>% mutate(sstmax = raster::extract(sstmax_lay,cell))
obiv <-  obiv %>% mutate(sstmean = raster::extract(sstmean_lay,cell))
obiv <-  obiv %>% mutate(sstrange= raster::extract(sstrange_lay,cell))

obiv <-  obiv %>% mutate(ss_salinity = raster::extract(ss_sal,cell))

obiv <-  obiv %>% mutate(bmtmin = raster::extract(btmin_min,cell))
obiv <-  obiv %>% mutate(bmtmax = raster::extract(btmax_min,cell))
obiv <-  obiv %>% mutate(bmtmean = raster::extract(bt_min,cell))
obiv <-  obiv %>% mutate(bmtrange= raster::extract(btr_min,cell))


# Filter out occurrences from localities with no environmental data
obiv <-  obiv %>% mutate(matchwithtemp = 1)
obiv$matchwithtemp[which(is.na(obiv$temperaturerange))] <- 0
# table(obiv$matchwithtemp)
# table(obiv$family[which(obiv$matchwithtemp==0)])

# remove non-matching, 231,308 remain:

obiv <-  subset(obiv,matchwithtemp == 1)

### Binary mineral columns

obiv <- obiv %>% mutate(isaragonite = 0)
obiv$isaragonite[which(obiv$mineral=="ara only")] <- 1
obiv <- obiv %>% mutate(iscalcite = 0)
obiv$iscalcite[which(obiv$mineral=="cal only")] <- 1
table(obiv$iscalcite)
obiv <- obiv %>% mutate(ismixed = 0)
obiv$ismixed[which(obiv$mineral=="mixed")] <- 1
table(obiv$ismixed)
table(obiv$isaragonite)

obiv$salinity <- as.numeric(obiv$salinity)
obiv$temperature <- as.numeric(obiv$temperature)
obiv$temperaturerange <- as.numeric(obiv$temperaturerange)

######### Filter based on depth and salinity
obiv <- obiv %>% filter(depth <= 100 & depth >= - 5)
obiv <- obiv %>% filter(BO_bathymax >= -200 )
# 121,254 occurrences
obiv <- obiv %>% filter(ss_salinity >= 30 & ss_salinity <= 40)
obiv <- obiv %>% filter(salinity >= 30 & salinity <= 40)
# 77,131
#
######################################################
#####
##### Add Biogeographic realms and grid cells
#####
#####
###
### GRID 8

library(icosa)
grid8 <- hexagrid(c(8,8), sp = F) # 
#plot(grid6)

# Locate cells
cells8 <- locate(grid8, cbind(obiv$decimalLongitude,obiv$decimalLatitude), randomborder = T)
obiv <- obiv %>% mutate(cells8 = cells8)
obiv$cells8 <- cells8

# unique cells (2950)
ocd8_unique <- unique(cells8)

n_cells8 <- sapply(ocd8_unique, function(x) length(which(obiv$cells8==x)))
#hist(n_cells8)

# get the mid longitudes and latitudes of the cells
longLat8 <- CarToPol(grid8@faceCenters, norad=T)
cellind8 <- sapply(ocd8_unique, function(x) which(rownames(longLat8)==x))
gridlon8 <- longLat8[cellind8,1]
gridlat8 <- longLat8[cellind8,2]
#plot(gridlon8,gridlat8,bg = rgb(0,0,0,0.9), pch = 21, col = NA, ylim = c(-90,90)) #, xlim = c(-180,180))

obiv <- obiv %>% mutate(gridlon8 = NA)
obiv$gridlon8 <- gridlon8[(sapply(obiv$cells8, function(x) which(ocd8_unique==x)))]
obiv <- obiv %>% mutate(gridlat8 = NA)
obiv$gridlat8 <- gridlat8[(sapply(obiv$cells8, function(x) which(ocd8_unique==x)))]


######################################################################################
###
### GRID 6


library(icosa)
grid6 <- hexagrid(c(6,6), sp = F) # 
#plot(grid6)

# Locate cells
cells6 <- locate(grid6, cbind(obiv$decimalLongitude,obiv$decimalLatitude), randomborder = T)
obiv <- obiv %>% mutate(cells6 = cells6)
obiv$cells6 <- cells6

# unique cells (2950)
ocd6_unique <- unique(cells6)

n_cells6 <- sapply(ocd6_unique, function(x) length(which(obiv$cells6==x)))
#hist(n_cells6)

# get the mid longitudes and latitudes of the cells
longLat6 <- CarToPol(grid6@faceCenters, norad=T)
cellind6 <- sapply(ocd6_unique, function(x) which(rownames(longLat6)==x))
gridlon6 <- longLat6[cellind6,1]
gridlat6 <- longLat6[cellind6,2]
#plot(gridlon6,gridlat6,bg = rgb(0,0,0,0.9), pch = 21, col = NA, ylim = c(-90,90)) #, xlim = c(-160,160))

obiv <- obiv %>% mutate(gridlon6 = NA)
obiv$gridlon6 <- gridlon6[(sapply(obiv$cells6, function(x) which(ocd6_unique==x)))]
obiv <- obiv %>% mutate(gridlat6 = NA)
obiv$gridlat6 <- gridlat6[(sapply(obiv$cells6, function(x) which(ocd6_unique==x)))]

######################################################################################
###
### GRID 4


library(icosa)
grid4 <- hexagrid(c(4,4), sp = F) # 
#plot(grid4)

# Locate cells
cells4 <- locate(grid4, cbind(obiv$decimalLongitude,obiv$decimalLatitude), randomborder = T)
obiv <- obiv %>% mutate(cells4 = cells4)
obiv$cells4 <- cells4

# unique cells (2950)
ocd4_unique <- unique(cells4)

n_cells4 <- sapply(ocd4_unique, function(x) length(which(obiv$cells4==x)))
#hist(n_cells4)

# get the mid longitudes and latitudes of the cells
longLat4 <- CarToPol(grid4@faceCenters, norad=T)
cellind4 <- sapply(ocd4_unique, function(x) which(rownames(longLat4)==x))
gridlon4 <- longLat4[cellind4,1]
gridlat4 <- longLat4[cellind4,2]
#plot(gridlon4,gridlat4,bg = rgb(0,0,0,0.9), pch = 21, col = NA, ylim = c(-90,90)) #, xlim = c(-160,160))

obiv <- obiv %>% mutate(gridlon4 = NA)
obiv$gridlon4 <- gridlon4[(sapply(obiv$cells4, function(x) which(ocd4_unique==x)))]
obiv <- obiv %>% mutate(gridlat4 = NA)
obiv$gridlat4 <- gridlat4[(sapply(obiv$cells4, function(x) which(ocd4_unique==x)))]


#### Assign biogeographic realms


require(sf)
require(rnaturalearth)
require(rgeos)

# meow
meow <- sf::read_sf("Meow/meow_ecos.shp") # it read as #c("MEOW.ECOREGION","MEOW.PROVINCE","MEOW.REALM","MEOW.Lat_Zone")

regions <- readOGR("Meow/meow_ecos.shp")
### Create province polygons
#Unite the spatial polygons for each region into one
# provinces <- unionSpatialPolygons(regions, regions$PROVINCE)
#Make a data frame that will have Province level info and above
# prov_data <- regions@data %>%
#  group_by(PROVINCE) %>%
#  summarise(PROV_CODE = PROV_CODE[1], REALM = REALM[1], RLM_CODE=RLM_CODE[1], Lat_Zone=Lat_Zone[1])

#merge the polygons with the new data file
#note the row.names argument to make sure they map to each other
# provinces <- SpatialPolygonsDataFrame(provinces,
#                                      data=data.frame(
#                                       join(data.frame(PROVINCE=names(provinces)),
#                                           prov_data),
#                                     row.names=row.names(provinces)))

## Joining by: PROVINCE

#######
#make realms shapefile
########
#make spatial polygons for realms
realms <- unionSpatialPolygons(regions, regions$REALM)

realmnames <- names(realms)

coords <- as.data.frame(obiv[,c("decimalLongitude","decimalLatitude")])
coordinates(coords) <- ~ decimalLongitude + decimalLatitude
proj4string(coords) <- proj4string(realms)

obiv <- obiv %>% mutate(realm = NA)
for(i in 1:12) {
  obiv$realm[which(!(is.na(sp::over(coords,realms[i]))))] <- i
}


# Assign entries in a grid cell southwest of Africa
obiv$realm[which(obiv$cells6 == "F9947")] <- 9

### Reorder the realm numbers into the sequence shown on the map

old <- 1:12
realms_reihenfolge <-  c(1,8,4,12,11,3,2,9,10,6,5,7)

oldrealms <- obiv$realm
oldrealms[oldrealms %in% old] <- realms_reihenfolge[match(oldrealms, old, nomatch = 0)]

table(oldrealms)

obiv$realm <- oldrealms

# check if the assignmentworked
table(obiv$realm)

# saveRDS("OBIS_bivalves_epifaunal_cleaned_2021_08_16_realms_and_grid.rds")
#######################################################################


