# Solution to the exercise: Brachiopods in OBIS 
# 2024-01-23
# https://adamtkocsis.com/rkheion/Exercises/2023-04-05_obis_brachiopods.html

# set working directory
setwd("/mnt/sky/Dropbox/Teaching/FAU/RConsolidation/2025/2025-11-10/obis_brachiopods/")

########################################----------------------------------------
# 1. Download this file: obis_brachiopoda.csv and put it in data

########################################----------------------------------------
# 2. read in the occurrence file (standard csv)
occs <- read.csv("data/obis_brachiopoda.csv")

########################################----------------------------------------
# 3. List out the genus entries
# genus entries
unique(occs$genus) # there are some that are ""

# This line returns the subset of the data where the genus column is not ""
filtered <- occs[which(!occs$genus==""), ]
# filtered <- occs[which(occs$genus!=""),] # literally the same result


########################################----------------------------------------
# 4. Continuing from where we left off: also omit rows where the species is ""
filterSP <- filtered[which(!filtered$species==""),]
unique(filterSP$species) # returns all the unique species names, the "" should not be in it!


########################################----------------------------------------
# 5. Species of Lingula genus

# subset of the Lingula occurrences
lingula <- filterSP[which(filterSP$genus=="Lingula"), ]
# subset(filterSP, filterSP$genus=="Lingula") # the same - it is somewhat more unstable

# counting the number of different species
# Any of these will do:
length(unique(lingula$species))
length(table(lingula$species))

# this is not good: it returns the number of rows in the lingula data.frame
# length(lingula$species)


########################################----------------------------------------
# 6. Plotting the coordinates
library(chronosphere)
# you might have to install the sf and chronosphere packages
# install.packages(c("sf","chronosphere"))

# download data from the internet
ne <- chronosphere::fetch("NaturalEarth")

# at any point you can create a new graphical device to plot into
# x11() # windows and linux
# quartz() # mac

# Basic plot of the landmassess (uses longitude/latitude coordinates)
plot(ne$geometry, col="gray")

# Put the points on the map.
# x coordinates are longitude,
# y coordiantes are latitude
points(
	x=lingula$decimallongitude,
	y=lingula$decimallatitude, col="red", pch="+"
)

# This is not good: (reversed coordinates)
## points(
## 	y=lingula$decimallongitude,
## 	x=lingula$decimallatitude, col="red", pch="+"
## )


 
########################################----------------------------------------
# 7. Number of species in every genus

# List of genus entries in the filterSP set (no "" species or genera!)
gen <- unique(filterSP$genus)

# set up a container to put the results in
# (one value for every genus)
nSpecies <- rep(NA, length(gen))
names(nSpecies) <- gen

# go through every genus (the index of the genus in the vector)
for(i in 1:length(gen)){
	# focal genus
	# focalGenus <- "Fosteria"
	focalGenus <- gen[i]
	
	# selecting data that pertain to focal genus
	speciesData <- filterSP[which(filterSP$genus==focalGenus), ]

	# count the number of species in the subset
	nSpecies[i] <- length(unique(speciesData$species))
}

########################################----------------------------------------
# 8. Make a histogram of the results
hist(nSpecies, breaks=6)

