## ATELIER DU CSBQ SUR L'UTILISATION DES LOGICIELS LIBRES POUR L'ANALYSE SPATIALE AVANCÉE

install.packages('rgdal','sp','spgrass6','RPostgreSQL','adehabitatHR','RQGIS')

library('rgdal')
library('sp')

## IMPORTATION D'UN FICHIER SHAPEFILE DANS R
#Choisissez le fichier meow_ecos.shp
meow=readOGR(file.choose(),layer='meow_ecos')
plot(meow)
spplot(meow,zcol='Lat_Zone')
sp.theme(set=TRUE,regions=list(col=terrain.colors(20)))
spplot(meow,zcol='Lat_Zone')

## IMPORTATION D'UN FICHIER RASTER DANS R
#Choisissez world_bedrock_topo.tif
world_topo=readGDAL(file.choose())
image(world_topo)

## CREATION D'UN SpatialPointsDataFrame géoréférencé
xy=rbind(c(-73.4,45.8),c(46.80,-19.2))
att=data.frame(villes=c('Montreal','Madagascar'))
mespoints=SpatialPointsDataFrame(coords=xy,data=att)
proj4string(mespoints)=CRS("+init=epsg:4326")
plot(meow)
plot(mespoints,add=TRUE,col=2)

## IMPORTATION D'UN FICHIER POSTGIS DANS R
library(rgdal)
dsn="PG:dbname='spatialdb' host='localhost' port='5432' user='postgres' password='yourpassword'"
ogrListLayers(dsn)
meow=readOGR(dsn,"meow")

## POUR DES REQUÊTES QUI RETOURNENT DES TABLES NON-SPATIALES
library(RPostgreSQL)
drv <- dbDriver("PostgreSQL")
con <- dbConnect(drv, user="postgres", password="yourpassword", dbname="spatialdb");
dbListTables(con)
# Espèces observées à l'intérieur de 1 degré des côtes du Canada
canada_whales <- dbGetQuery(con,"SELECT a.* FROM obis_whales a, tm_world b WHERE ST_DWithin(a.geom,b.geom,1) AND b.name='Canada'");
unique(canada_whales$tname)

## COMBINAISON DES DEUX

dbGetSp <- function(dbInfo,query) {
  if(!require('rgdal')|!require(RPostgreSQL))stop('missing rgdal or RPostgreSQL')
  d <- dbInfo
  tmpTbl <- sprintf('tmp_table_%s',round(runif(1)*1e5))
  dsn <- sprintf("PG:dbname='%s' host='%s' port='%s' user='%s' password='%s'",
                 d$dbname,d$host,d$port,d$user,d$password
  )
  drv <- dbDriver("PostgreSQL")
  con <- dbConnect(drv, dbname=d$dbname, host=d$host, port=d$port,user=d$user, password=d$password)
  tryCatch({
    sql <- sprintf("CREATE UNLOGGED TABLE %s AS %s",tmpTbl,query)
    res <- dbSendQuery(con,sql)
    nr <- dbGetInfo(res)$rowsAffected
    if(nr<1){
      warning('There is no feature returned.');
      return()
    }
    sql <- sprintf("SELECT f_geometry_column from geometry_columns WHERE f_table_name='%s'",tmpTbl)
    geo <- dbGetQuery(con,sql)
    if(length(geo)>1){
      tname <- sprintf("%s(%s)",tmpTbl,geo$f_geometry_column[1])
    }else{
      tname <- tmpTbl;
    }
    out <- readOGR(dsn,tname)
    return(out)
  },finally={
    sql <- sprintf("DROP TABLE %s",tmpTbl)
    dbSendQuery(con,sql)
    dbClearResult(dbListResults(con)[[1]])
    dbDisconnect(con)
  })
}
coninfo=list(host='localhost',dbname='spatialdb',port='5432',user='postgres',password='yourpassword')
#Sélectionner toutes les observations dans le Golfe du Mexique
vec<-dbGetSp(coninfo,"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'")

plot(vec)

## Calcul des densités de points (kernel density) à partir des résultats de PostGIS
library(adehabitatHR)
library(raster)
# Densité des observations dans le golfe du Mexique.
ku=kernelUD(vec,grid=1000,extent=0.3,h=0.5)
gv=getverticeshr(ku)
image(ku) # Couche raster
plot(vec,pch=1,cex=0.1,alpha=0.1,add=TRUE)
plot(gv,add=TRUE)
out=raster(ku)

## UTILISATION DES FONCTIONS QGIS DANS R
library(RQGIS)
my_env <- set_env()
dir_tmp <- '/home/glaroc/Workshops/specialspatial/' #Remplacer par un répertoire existant
find_algorithms(search_term = "polygon", 
                qgis_env = my_env)
params <- get_args_man(alg = "qgis:polygoncentroids", qgis_env = my_env) #Retrouver la liste des paramètres
params$INPUT_LAYER  <- meow
# destination du fichier shapefile en sortie
params$OUTPUT_LAYER <- file.path(dir_tmp, "ger_coords.shp")
out <- run_qgis(alg = "qgis:polygoncentroids",
                params = params,
                load_output = params$OUTPUT_LAYER,
                qgis_env = my_env)

## UTILISATION DES FONCTIONS GRASS DANS R
library(spgrass6)
## Linux
#initGRASS(gisBase='/usr/lib/grass70',gisDbase='/home/glaroc/OSGIS/GRASS',location='monde',mapset='baleines',override=TRUE)
## Windows
## Remplacer gisBase et gisDbase par les dossiers appropriés. 
initGRASS(gisBase='C:/OSGeo4W64/apps/grass/grass-7.0/',gisDbase='/Users/yourname/GIS/GRASS',location='monde',mapset='baleines',override=TRUE)
## Mac
## initGRASS(gisBase='/usr/lib/grass70',gisDbase='/home/yourname/GIS/GRASS',location='monde',mapset='baleines',override=TRUE)

gmeta6()
execGRASS("g.list",type="vect") # Liste des fichiers vecteur dans GRASS
execGRASS("g.list",type="rast") # Liste des fichiers raster
pays<-readVECT6('pays')
out=execGRASS("v.extract",flags=c("overwrite"),input="pays", output="Canada",where="\"NAME='Canada'\"")
canada<-readVECT6('Canada')
plot(canada)
execGRASS("g.region",vect="Canada",res="0.1")
execGRASS("v.to.rast",input="Canada",output="Canada_rast",use="val",flags=c("overwrite"))
execGRASS("r.mapcalculator",amap="world_topo",bmap="Canada_rast",formula="A*B",outfile="topo_canada")
execGRASS("g.remove",rast="topo_canada")
topo_canada<-readRAST6('topo_canada')
image(topo_canada)
hist(topo_canada@data[,1])