---
title: "Moss map occurrences"
output: html_notebook
editor_options: 
  chunk_output_type: console
---

```{r}

library(raster)
library(tidyverse)
library(data.table)
library(sf)
library(lubridate)
```

```{r}

chelsa_bio<-raster("/Users/tdg139/Library/CloudStorage/Dropbox/Postdoc/other manuscripts/moss_revision/data/CHELSA_bio1_1981-2010_V.2.1.tif")

hypnum<-fread("data/hypnum.csv")
hypnum<-hypnum %>% select(species,decimalLatitude,decimalLongitude,coordinateUncertaintyInMeters,basisOfRecord)

sanionia<-fread("data/sanionia.csv")
sanionia<-sanionia %>% select(species,decimalLatitude,decimalLongitude,coordinateUncertaintyInMeters,basisOfRecord)

both<-full_join(hypnum, sanionia)
both<-both %>% filter(basisOfRecord %in% c("HUMAN_OBSERVATION", "OCCURRENCE", "OBSERVATION")) %>% drop_na(decimalLatitude,decimalLongitude) %>% filter(coordinateUncertaintyInMeters<=1000)

both_sp<-st_as_sf(both, coords=c("decimalLongitude","decimalLatitude"),crs=4326,remove = F)

#extract
mossclim<-data.frame(both_sp, "antemp"=raster::extract(chelsa_bio,both_sp,cellnumbers=TRUE))

mossclim<-mossclim %>% group_by(species) %>% distinct(antemp.cells,.keep_all=T)

mossclim2<-mossclim %>% mutate(antemp=antemp.CHELSA_bio1_1981.2010_V.2.1*0.1+-273.15)

mossclim2<-st_as_sf(mossclim2, coords=c("decimalLongitude","decimalLatitude"),crs=4326,remove = F)

ggplot(mossclim2)+geom_sf(europe)+geom_sf(aes(col=species))

#Species means in annual temperature
mossclim2 %>% group_by(species) %>% summarise(mean(antemp,na.rm=T))

#Occurrence map
library(rnaturalearthdata)
library(rnaturalearth)
world2<-ne_download(scale=110,type="countries")
europe <- world2[world2$REGION_UN=="Europe",]
europe<-europe %>% filter(!SOVEREIGNT=="Russia")

# make a fishnet grid over the countries
grd <- st_make_grid(europe, cellsize = .5)
# find which grid points intersect `polygons` (countries) 
# and create an index to subset from
index <- which(lengths(st_intersects(grd, europe)) > 0)
# subset the grid to make a fishnet
fishnet <- grd[index]

map1<-fishnet |>  
  st_as_sf() |> # cast to sf
  mutate(grid_id = row_number()) |> # create unique ID
  st_join(mossclim2) %>% distinct(grid_id,species,.keep_all=T)
 
 map1<-map1 %>% group_by(grid_id) %>% mutate(counts=n())
 map1<-map1 %>% mutate(species2=ifelse(counts==2, "Both",species))

filter_sf <- function(.data, xmin = NULL, xmax = NULL, ymin = NULL, ymax = NULL) {
  bb <- sf::st_bbox(.data)
  if (!is.null(xmin)) bb["xmin"] <- xmin
  if (!is.null(xmax)) bb["xmax"] <- xmax
  if (!is.null(ymin)) bb["ymin"] <- ymin
  if (!is.null(ymax)) bb["ymax"] <- ymax
  sf::st_filter(.data, sf::st_as_sfc(bb), .predicate = sf::st_within)
}

map2<-map1 %>%
  filter_sf(ymin = 25) %>% 
ggplot()+geom_sf(aes(fill=species2))+ scale_fill_manual(values=c("darkseagreen","darkgoldenrod2","darkorange3","grey"))+labs(fill="Species")+theme_bw()

ggsave("mossspecies_map_lowerres.png",width=6,height=5)

```

