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

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

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

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

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

RPostGIS <- function(coninfo,query) {
  dsn=paste("PG:dbname='",coninfo$dbname,"' host='",coninfo$host,"' port='",coninfo$port,"' user='",coninfo$user,"' password='",coninfo$password,"'", sep='')
  drv <- dbDriver("PostgreSQL")
  con <- dbConnect(drv, user=coninfo$user, password=coninfo$password, dbname=coninfo$dbname)
  res <- dbSendQuery(con,paste('CREATE TABLE tmp1209341251dva1 AS ',query,sep=''))
  geo <- dbGetQuery(con,"SELECT f_geometry_column from geometry_columns WHERE f_table_name='tmp1209341251dva1'");
  if(length(geo)>1){
    tname=paste("tmp1209341251dva1(",geo$f_geometry_column[1],")")
  }else{
    tname="tmp1209341251dva1"
  }
  out <- tryCatch(readOGR(dsn,tname), finally=dbSendQuery(con,'DROP TABLE tmp1209341251dva1'))
  dbDisconnect(con)
  return(out)
}

coninfo=list(host='localhost',dbname='spatialdb',port='5432',user='postgres',password='yourpassword')
vec<-RPostGIS(coninfo,"SELECT a.id_0, a.geom_pt, a.tname, a.datecollec FROM obis_whales a, meow_ecos 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)
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)
plot(pol,add=TRUE)
out=raster(ku)

## UTILISATION DES FONCTIONS GRASS DANS R
library(spgrass6)
## Linux
#initGRASS(gisBase='/usr/lib/grass64',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-6.4.3/',gisDbase='/Users/yourname/GIS/GRASS',location='monde',mapset='baleines',override=TRUE)
## Mac
## initGRASS(gisBase='/usr/lib/grass64',gisDbase='/home/yourname/GIS/GRASS',location='monde',mapset='baleines',override=TRUE)

gmeta6()
execGRASS("g.list",type="vect")
execGRASS("g.list",type="rast")
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])