[R] Simplifier un fond de carte avec R et PostGIS

Rpostgisgeneralisation

Ce billet fait suite au précédent billet intitulé « Combiner R et PostGIS, quelle bonne idée !« . Ce billet visait à montrer comment utiliser conjointement les points forts de ces deux outils dans un environnement de développement cohérent [VOIR].  Dans cette nouvelle publication, je montre comment simplifier un fond de carte en m’appuyant sur cette architecture (R+PostGIS). Plus précisément,  je montre comment on peut travailler avec la puissance des fonctionnalités de PostGIS tout en restant dans l’environnement de R. Pour cela, je mets en oeuvre la méthode de généralisation cartographique présentée sur le site http://trac.osgeo.org/ [VOIR]

Avant de commencer

Pour commencer, il faut que R et PostGIS soient installés sur votre système. Vous devez aussi télécharger le script ci-dessous permettant d’échanger les données entre R et PostGIS.

download

Le fond de carte utilisé dans l’exemple ci-dessous est téléchargeable sur ce blog dans la rubrique ressources : [VOIR]

C’est parti

##########################################################################
# GENERALISATION CARTOGRAPHIQUE AVEC POSTGIS ET R
# d'après http://trac.osgeo.org/postgis/wiki/UsersWikiSimplifyPreserveTopology
# Nicolas LAMBERT, janvier 2015
##########################################################################

# [STEP 0] Configuration et appel des librairies

setwd("votre-repertoire-de-travail") # ici votre répertoire de travail
source("sources/Rpostgis.R")
con <- dbConnect(dbDriver("PostgreSQL"), dbname="test") # ici vos paramètres de connexion

# [STEP 1] Import du fond de carte dans PostGIS

# Shapefile -> spatial dataframe
fdc <- readOGR(dsn ="geom",layer = "GEOM_EU34_NUTS2010")
countries <-fdc[fdc@data$level0==1,]
# On met les noms des champs en minuscule pour eviter tout problème dans PostGIS
names(countries) <- tolower(names(countries))
# Spatial dataframe -> PostGIS
dbWriteSpatial(con, countries, schemaname="public", tablename="countries", replace=T,srid=3035)

# [STEP 2] Définition des paramètres de généralisation

# Nom de la couche en entrée
input_layer <- "countries"
# Chanp unique identifiant les unités territoriales
id <- "id"
# Niveau de simplification
simplification <- 100000
# Nom de la couche de sortie
output_layer <- "countries_generalized"

# [STEP 3] Multipolygons -> polygons en gardant les identifiants

# Eriture de la requète PostGIS
q<-paste("create table poly as (select ",id,", (st_dump(geom)).*from ",input_layer,");",sep="")
# Execution de la requete
res<-dbSendQuery(con,q)
# Libération des ressources
dbClearResult(res)

# [STEP 4] Extractions des contours exterieurs des polygones (rings)

q<-paste("create table rings as (select ",id,", st_exteriorRing((st_dumpRings(geom)).geom) as g from poly);",sep="")
res<-dbSendQuery(con,q)
dbClearResult(res)

# [STEP 5] Union de toutes les lignes et simplification des contours (paramétrable)

q<-paste("create table simplerings as (select st_simplifyPreserveTopology(st_linemerge(st_union(g)), ",simplification,") as g from rings);",sep="")
res<-dbSendQuery(con,q)
dbClearResult(res)

# [STEP 6] Extraction des lignes en tant qu'objets individuels (multi-lines -> lines)

q<-"create table simplelines as (select (st_dump(g)).geom as g from simplerings);"
res<-dbSendQuery(con,q)
dbClearResult(res)

# [STEP 7] Reconstruction des polygons à partir des lignes simplifiées

q<-"create table simplepolys as (select (st_dump(st_polygonize(distinct g))).geom as g from simplelines);"
res<-dbSendQuery(con,q)
dbClearResult(res)

# [STEP 8] Jointure spatiale pour récupérer les id perdus pendant la bataille

# Ajout de l'id
q<-"alter table simplepolys add column gid serial primary key;"
res<-dbSendQuery(con,q)
dbClearResult(res)
# création de l'index
q<-"create index simplepolys_geom_gist on simplepolys  using gist(g);"
res<-dbSendQuery(con,q)
dbClearResult(res)
# Recuperation des id par jointure spatiale
q<-paste("create table simple_tmp as (select ",id,", g from ",input_layer," d, simplepolys s  where st_contains(d.geom, st_pointOnSurface(s.g)));",sep="")
res<-dbSendQuery(con,q)
dbClearResult(res)

# [STEP 9] Recréation du fond de carte

q<-paste("create table ",output_layer," as (select ",id,", st_collect(g) as geom  from simple_tmp group by ",id," );",sep="")
res<-dbSendQuery(con,q)
dbClearResult(res)

# [STEP 10] Import et affichage dans R

# PotsGIS -> R
output <- dbReadSpatial(con, schemaname="public", tablename=output_layer, geomcol="geom")
# Affichage
par(mfrow=c(1,2))
plot(countries,col="#AAAAAA")
title(main="AVANT")
plot(output,col="#AAAAAA")
title(main="APRES")
par(mfrow=c(1,1))

# [STEP 11] Nettoyage de toutes les tables temporaires créées dans la base PostGIS

res<-dbSendQuery(con,"SELECT DropGeometryTable('rings')")
dbClearResult(res)
res<-dbSendQuery(con,"SELECT DropGeometryTable('poly')")
dbClearResult(res)
res<-dbSendQuery(con,"SELECT DropGeometryTable('simple_tmp')")
dbClearResult(res)
res<-dbSendQuery(con,paste("SELECT DropGeometryTable('",output_layer,"')",sep=""))
dbClearResult(res)
res<-dbSendQuery(con,"SELECT DropGeometryTable('simplelines')")
dbClearResult(res)
res<-dbSendQuery(con,"SELECT DropGeometryTable('simplepolys')")
dbClearResult(res)
res<-dbSendQuery(con,"SELECT DropGeometryTable('simplerings')")
dbClearResult(res)
res<-dbSendQuery(con,"SELECT DropGeometryTable('countries')")
dbClearResult(res)

# [STEP 12] Si besoin, conversion en shapefile

writeOGR(output, ".", "output", driver="ESRI Shapefile")

Résultats

result


Vous aimerez aussi...

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *