[R] Combiner R et PostGIS, quelle bonne idée !
Pourquoi Combiner R et PostGIS ?
Combiner R et PostGIS permet d’utiliser dans un environnement de travail intégré regroupant à la fois la puissance de R (gestion des données, traitements statistiques, …) et la puissance de PostGIS (géotraitements, requêtes spatiales, …) sans rupture de charge. Cette alliance permet notamment de : traiter, analyser et cartographier les données (R) ; bénéficier des dernières évolutions de PostGIS 2 (gestion topologique, indexation, géotraitements, …) ; structurer les données de façon efficace (PostGIS) ; favoriser le recherche reproductible (R & PostGIS). Et tout cela avec les avantages du libre et du gratuit (nombreuses mises à jour, grosse communauté d’utilisateurs). Bref, Cette combinaison permet d’offrir un environnement de développement complet et puissant adapté à la très grande majorité des traitements utiles en géographie, en statistique et géomatique. Que demander de plus…
Prérequis
Pour commencer, il faut faire quelques installations :
installer R (http://www.r-project.org/)
installer et configurer PostGIS (http://postgis.net/install/ & http://postgis.net/docs/postgis_installation.html)
Comment transférer les données de R à PostGIS et inversement ?
Théoriquement, il est possible d’importer et d’exporter des spatialdataframes dans PostGis avec le package rgdal. Néanmoins, après de nombreux essais infructueux, une solution basée sur l’import/export des géométries au format WKT s’avère nettement plus efficace. La meilleure solution trouvée sur le web est proposée sur le site de Philipp Hunziker. Sur ce site, l’auteur propose deux fonctions permettant d’importer ou d’exporter les géométries et les tables attributaires. Les fonctions proposées sur ce site ne gérant ni les projections ni l’import de polygones uniques, je propose ci-dessous une version un petit peu mise mise à jour.
Une fois le script Rpostgis.R téléchargé et sourcé, la fonction dbWriteSpatial() permet de transférer un spatialdataframe vers PostGIS.
La fonction dbReadSpatial() permet de de transférer un une table PostGIS vers un spatialdataframe
Exemple d’utilisation :
# Paramètres de connexion con <- dbConnect(dbDriver("PostgreSQL"), dbname="test") # Shapfile -> R -> PostGIS df <- readOGR(dsn ="geom/",layer = "polys") dbWriteSpatial(con, df, schemaname="public", tablename="newtable", replace=T,srid=3035) # PostGIS -> R spdf <- dbReadSpatial(con, schemaname="public", tablename="newtable", geomcol="geom") proj4string(spdf) head(spdf@data) plot(spdf)
Exécuter une requête PostGIS dans R
L’exécution de requêtes PostGIS dans R s’effectue avec le package RPostgreSQL
# Appel du package RPostgreSQL require(RPostgreSQL) # Définition des paramètres de connexion con <- dbConnect(dbDriver("PostgreSQL"), dbname="test") # Requête (e.g zone tampon) q<-"CREATE TABLE nom-de-ma-nouvelle-table-postgis- AS SELECT id,ST_Buffer(geom,200000) AS geom FROM nom-de-ma-table-postgis" # Envoi de la requête vers PostGIS res<-dbSendQuery(con,q) # Libération des ressources pour effectuer éventuellement une autre requête dbClearResult(res)
Un exemple de traitement simple (e.g. Polygons -> polylines)
source("sources/Rpostgis.R") # 1. Connexion con <- dbConnect(dbDriver("PostgreSQL"), dbname="test") # 2. Transfert d'une couche sur les pays européens dans PostGIS df <- readOGR(dsn ="data/",layer = "NUTS0_RG_20M_2006") dbWriteSpatial(con, df, schemaname="public", tablename="pays", replace=T,srid=3035) # 3. Exécution d'un requête (Polygones -> polylignes) q<-"CREATE TABLE bound AS SELECT id,ST_Boundary(geom) AS geom FROM pays" res<-dbSendQuery(con,q) dbClearResult(res) # 4. Récupération du résultat dans R df2 <- dbReadSpatial(con, schemaname="public", tablename="bound", geomcol="geom") # 5. Affichage par(mfrow=c(1,2)) plot(df,col="red") plot(df2,col="red") par(mfrow=c(1,1)) # 6. Nettoyage de la base PostGIS res1<-dbSendQuery(con,"SELECT DropGeometryTable('bound')") dbClearResult(res1) res2<-dbSendQuery(con,"SELECT DropGeometryTable('pays')") dbClearResult(res2)
Conclusion
Si l’exemple montré ici est simple (il aurait pu être effectué intégralement dans R sans PostGIS), il montre néanmoins un aperçu des possibilités offertes par l’articulation R<->PostGIS ainsi que la logique d’une telle utilisation. Dans un prochain billet, je montrerai comment utiliser cet environnement de développement pour réaliser une généralisation cartographique qui conserve la topologie.
Bonjour,
Merci beaucoup pour ce code.
J’ai toutefois du faire des modifications aussi bien dans dbWriteSpatial que dans dbReadSpatial pour que ces fonctions fonctionnent.
Dans dbWriteSpatial, j’ai modifié la ligne 43 : query5 <- paste("SELECT UpdateGeometrySRID('",tablename,"','geom',",srid,")",sep="")
par query5 <- paste("SELECT UpdateGeometrySRID('",schemaname,"','",tablename,"','geom',",srid,")",sep="")
Dans dbReadSpatial j'ai modifié schema.table = paste(schemaname, ".", tablename, sep="") en schema.table = paste0("\"",schemaname, "\".\"", tablename,"\"") ainsi que query1 <- paste("SELECT ST_SRID(geom) FROM",tablename,"LIMIT 1")
en query1 <- paste("SELECT ST_SRID(\"",tablename,"\".geom) FROM \"",schemaname,"\".\"",tablename,"\" LIMIT 1",sep="")
Cordialement
Serge