## QCBS WORKSHOP ON THE USE OF OPEN SOURCE TOOLS FOR ADVANCED SPATIAL ANALYSIS
## CREATOR: Guillaume Larocque (glaroc@gmail.com)

install.packages(c('rgdal','sp','rgrass7','RPostgreSQL','adehabitatHR','RQGIS3','rpostgis','sf', 'ggmap'))
library(rgdal)

## IMPORTING A SHAPEFILE IN R
# Choose the file meow_ecos.shp
meow=readOGR(file.choose(),layer='meow_ecos')
plot(meow)
spplot(meow,zcol='Lat_Zone')
spplot(meow,zcol='Lat_Zone',col.regions=terrain.colors(10))

## IMPORTING A RASTER FILE IN R
# Choose world_bedrock_topo. Warning, this may be slow!
world_topo=readGDAL(file.choose())
image(world_topo)

## CREATION OF A GEOREREFENCED SpatialPointsDataFrame
xy=rbind(c(-73.4,45.8),c(46.80,-19.2))
att=data.frame(cities=c('Montreal','Madagascar'))
mypoints=SpatialPointsDataFrame(coords=xy,data=att)
proj4string(mypoints)=CRS("+init=epsg:4326")
plot(meow)
plot(mypoints,pch=16,col="red",add=TRUE)

## IMPORTING A POSTGIS FILE IN R
library('rpostgis')
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv,host = "localhost", user="postgres", port="5432", password="your_password", dbname="spatialdb");
vec<-pgGetGeom(con, query = "SELECT a.id, a.geom, a.tname, a.datecollec FROM obis_whales a, meow b WHERE ST_Within(a.geom,b.geom) AND b.ecoregion='Northern Gulf of Mexico'")

## Calculate point densities (kernel density) from PostGIS results in previous step using AdehabitatHR package
library('adehabitatHR')
ku=kernelUD(vec,grid=1000,extent=0.3,h=0.5)
gv=getverticeshr(ku)
image(ku) # Raster layer
plot(vec,pch=1,cex=0.1,alpha=0.1,add=TRUE)
plot(gv,add=TRUE)
library(raster)
out=raster(ku)

## SF PACKAGE
library(sf)
#Convert between sp and sf
#mypoints_sf<-st_as_sf(mypoints)
#Read from PostGIS
meow<-st_read(con, "meow")
plot(meow["ecoregion"])
obis_whales<-st_read(con,"obis_whales")
# Points within Cold Temperate Northwest Atlantic ecoprovince
CTN<-st_read(con, query = "SELECT a.id, a.geom, a.tname, a.datecollec FROM obis_whales a, meow b WHERE ST_Within(a.geom,b.geom) AND b.province='Cold Temperate Northwest Atlantic' AND a.tname='Orcinus orca'")
plot(CTN)

# This is equivalent to:
A<-obis_whales[obis_whales$tname=='Orcinus orca',]
B<-meow[meow$province=='Cold Temperate Northwest Atlantic',]
vec2<-st_within(A,B,sparse=FALSE) #returns a logical index
ll<-rowSums(vec2) > 0L
CTN2<-A[as.matrix(ll),]
plot(CTN2['tname'])

# Add buffer around points
buf<-st_buffer(CTN2,dist=1) #1 degree buffer
plot(buf['tname'])
# Merge all buffers into a single polygon
bufu<-st_union(buf)
plot(CTN2['tname'], reset=FALSE)
plot(bufu, add=TRUE)

bb<-st_bbox(CTN2)
library(ggmap)
names(bb)<-c("left","bottom","right","top")
#Show points on a Open Street Maps background
map<-get_map(bb, zoom = 4, source="stamen") 
ggmap(map)+geom_point(data=as.data.frame(CTN2),aes(longitude,latitude))



#
# EXERCISE
#
# How many species of toothed whales are found within 50km of Ponta Delgada, 
# in the Azores, Portugal (lat:37.761672, long:-25.627855)? 
# Use the sf package to find the solution. You can use the st_transform() 
# function to convert to and from UTM Zone 26 N (EPSG 32626) 
# and lat/long (EPSG 4326)


#
# USING QGIS FUNCTIONS IN R
#
library('devtools')
devtools::install_github("jannes-m/RQGIS3")

library(RQGIS3)
my_env <- set_env()
dir_tmp <- '/home/glaroc/Workshops/specialspatial/' #Replace with an existing folder
library("RQGIS")
find_algorithms(search_term = "polygon", qgis_env = my_env)
params <- get_args_man(alg = "qgis:polygoncentroids", qgis_env = my_env) #Retrieve list of parameters
params$INPUT_LAYER  <- meow
# destination of the output shapefile 
params$OUTPUT_LAYER <- file.path(dir_tmp, "ger_coords.shp")
# Find centroids of polygons of marine ecoregions
out <- run_qgis(alg = "qgis:polygoncentroids",
                params = params,
                load_output = params$OUTPUT_LAYER,
                qgis_env = my_env)
plot(out)





## USING GRASS FUNCTIONS IN R
library(rgrass7)
## Linux
# initGRASS(gisBase='/usr/lib/grass74',gisDbase='/home/yourname/GIS/GRASS',location='newLocation',mapset='Ontario',override=TRUE)
## Windows
## Remplace gisBase and gisDbase by the appropriate folders. 
initGRASS(gisBase='C:/Program Files/QGIS 3.6/apps/grass/grass76',gisDbase='/Users/yourname/GIS/GRASS',location='world',mapset='whales',override=TRUE)
 #initGRASS(gisBase='C:/OSGeo4W64/apps/grass/grass-7.4.3/',gisDbase='/Users/yourname/GIS/GRASS',location='world',mapset='whales',override=TRUE)
## Mac
## initGRASS(gisBase='/usr/lib/grass64',gisDbase='/home/yourname/GIS/GRASS',location='newLocation',mapset='Ontario',override=TRUE)

#Files can be imported within the GRASS interface, using writeVEC or writeRAST or within QGIS with the GRASS plugin
gmeta()
execGRASS("g.list",type="vect")
execGRASS("g.list",type="rast")
tm_world<-readVECT('tm_world')
#Choose Canada from polygon map of countries
out=execGRASS("v.extract",flags=c("overwrite"),input="tm_world", output="Canada",where="\"NAME='Canada'\"")
#Import into R
canada<-readVECT('Canada')
plot(canada)
#Define the GRASS region as all of Canada with a resolution of 0.1 degrees
execGRASS("g.region",vector="Canada",res="0.1")
#Convert to raster
execGRASS("v.to.rast",input="Canada",output="Canada_rast",use="val",flags=c("overwrite"))
#Extract topographic map within Canada
execGRASS("r.mapcalc",expression="topo_canada=world_topo*Canada_rast")
execGRASS("g.region",vector="Canada",res="0.1")
#Read topo within Canada as R Raster object
topo_canada<-readRAST('topo_canada')
image(topo_canada)
hist(topo_canada@data[,1])