Species distribution modelling

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")

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

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

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

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

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

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)

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

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

## This is MaxEnt version 3.4.1

dist_now <- predict(me, worldclim) 

## This is MaxEnt version 3.4.1

names(worldclim_future) <- names(worldclim)
dist_future <- predict(me, worldclim_future) 

## This is MaxEnt version 3.4.1

dist_now[dist_now<0.3]<-NA
dist_future[dist_future<0.3]<-NA

library(rnaturalearth)
europe <- ne_countries(continent = 'europe',10)


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

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

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)

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