---
title: "Species distribution modelling"
author: "Duncan Golicher"
date: "01/12/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE,message = FALSE)
```

```{r}
library(dismo)
library(sf)

library(mapview)
library(tmap)
library(RPostgreSQL)
library(dplyr)
library(raster)
library(rgdal)
data(World)
mapviewOptions(basemaps =c("OpenStreetMap",
"Stamen.TonerLite", "Esri.WorldImagery", 
"Esri.WorldShadedRelief", 
"Stamen.Terrain"))

conn <- dbConnect("PostgreSQL", host = "postgis",
  dbname = "gis_course" ,user = "gis_course", 
  password = "gis_course123")
```



```{r}
query<-"select geometry,species from gbif where country = 'United Kingdom'"
records<-st_read(conn, query=query)
```


```{r}
query<-"select geometry,species from 
gbif where country = 'United Kingdom' 
order by random() limit 2000"

st_read(conn, query=query) %>% mapview(zcol="species", burst=TRUE) -> rec_map
rec_map
```

```{r}
gridded<-st_read(conn,query="select * from gbif_gridded")
gridded$occurence<-(gridded$count>0)*1 # Make a column of ones and zeros
mapview(gridded,zcol="species",burst=TRUE) -> grid_map
grid_map
```


```{r}
occurence <- st_read(conn,query="select * from gbif 
where genus = 'Hippocrepis' 
order by random() limit 10000")
```



```{r}
limits<-st_convex_hull(occurence)# Make a hull around the points 
limits<-st_buffer(limits,1) # Set the area of interest to 1 degree around the points

limits<-as(limits,"Spatial") # Needed as the code uses the older form of R spatial objects.
occurence <- as(occurence, 'Spatial') ## Change to the old spatial class
```

```{r}
path<- "~/geoserver/data_dir/rasters/worldclim"
res<-5
worldclim<-getData(name = "worldclim", var = "bio", res = 5,path=path) ## Get all the bioclim layers
future <- getData('CMIP5', var='bio', res=5, rcp=85, model='AC', year=50, path=path)
worldclim<-raster::crop(worldclim,limits)
worldclim_future<-raster::crop(future,limits)

```


```{r}
fold <- kfold(occurence, k=5)
occtest <- occurence[fold == 1, ]
occtrain <- occurence[fold != 1, ]

me <- maxent(worldclim, occtrain, path="maxent")
```

```{r}
dist_now <- predict(me, worldclim) 

names(worldclim_future) <- names(worldclim)
dist_future <- predict(me, worldclim_future) 
dist_now[dist_now<0.3]<-NA
dist_future[dist_future<0.3]<-NA
```


```{r}
library(rnaturalearth)
europe <- ne_countries(continent = 'europe',10)


qtm(dist_now) + 
  qtm(europe, fill=NULL) + 
  tm_legend(legend.outside=TRUE)
```


```{r}
qtm(dist_future) + 
  qtm(europe, fill=NULL) + 
  tm_legend(legend.outside=TRUE)
```


```{r}
query<-"select st_transform(geom,4326) geom 
from bedrock where rcs_d ilike 'CHALK%' 
or rcs_d ilike 'LIMESTONE%'"

chalk<- st_zm(st_read(conn, query=query))
bounds<-as(chalk,"Spatial")



dist_now<-raster::crop(dist_now, bounds)
dist_future<-raster::crop(dist_future, bounds)
qtm(dist_now) + qtm(europe, fill = NULL) + tm_legend(legend.outside=TRUE)
```


```{r}
d<-raster::mask(dist_now,chalk)
tmap_options(check.and.fix = TRUE) 
mapview(d)
```

