Export IFC vers PostGIS 2ème Partie

IFC2PG : Industry Foundation Classes vers PostGIS 2ème Partie

 

Introduction

Ce billet est la suite de la première partie que vous trouverez ici IFC vers PostGIS 1ère partie, prenez en connaissance avant de poursuivre votre lecture.

Le développement de l’outil est réalise dorénavant sous Visual Studio 2010 en VB.Net.

Nouveau composant

J’ai ajouté un nouveau composant DotSpatial qui permet de visualiser l’export 2D des objets du fichier IFC.

DotSpatial est développé par les membres de l'équipe open source SIG de MapWindow, et la communauté de développeurs d'OSGeo .NET

Evolution

J’ai rajouté les objets suivants dans l’export :

• IfcSpace
• IfcSlab
• IfcColumn
• IfcDoor
• IfcWindow

Ces objets sont exportés sous forme de géométrie de type polygon.

J’ai modifié l’export de l’objet IfcWall qui maintenant donne une géométrie de type polygon.

L’interface s’est enrichi d’une vue 2D permettant de visualiser le résultat des exports.

Installation

Téléchargez l’outil IFC2PG : ifc2pg.rar

Décompressez l’archive sous c:\

Créez un raccourci sur le bureau.

 

Utilisation

 

Lancez l’application par un double-clic sur le raccourci.

1. Arborescence du fichier IFC,
2. Zone d’affichage 3D du fichier IFC (IFCEngine OCX) ,
3. Gestion des niveaux (DoSpatial),
4. Zone d’affichage 2D des objets exportés (DotSpatial)

Cet icône permet d’ouvrir un fichier IFC.

Cet icône permet de se connecter à une base de données PostgreSQL/PostGIS.

Ouvrez le fichier qui se trouve dans le répertoire IFC

Export des données

Pour réaliser l’export des données, sélectionner un objet IfcBuildingStorey puis faites un clic droit pour faire apparaitre un menu contextuel.

Ce menu permet de réaliser un export par type d’objet (IfcWallStandardCase, IfcSpace, IfcSlab, IfcColumn, IfcDoor, IfcWindow) en sélectionnant la ligne correspondante dans le menu, ou bien un export par IfcBuildindStorey, tous les objets seront alors exportés.

Ce menu permet également de gérer l’affichage des objets dans la vue 3D par type d’objet ou pour tous les objets contenus par l’objet IfcBuildingStorey.

Choisissez Export de tous les Objets dans le menu, vous devriez obtenir ceci :

Pour que les objets soient transférés dans PostGIS, il faut se connecter à une base de données:

1. Cliquez sur le bouton OK,
2. Sélectionnez un objet IfcBuildingStorey,
3. faites un clic droit, puis choisissez Export de tous les Objets dans le menu,

Il ne vous reste plus qu’a allez visualiser le contenu tables avec par exemple OpenJump.

Important

Je ne vide pas les tables avant export, donc pensez à les vider pour ne pas vous retrouvez avec des doublons, je vais réfléchir à comment gérer ce risque de doublons.

Conclusion

Nous voila à la fin de cette seconde partie, je vais continuer à améliorer l’interface de l’outil, ajouter de nouveaux objets… Rendez-vous prochainement dans une troisième partie.

Je vais également réfléchir à un export 3D sous PostGIS 2.0

Export IFC vers PostGIS 1ere Partie

 

IFC2PG : Industry Foundation Classes vers PostGIS 1ere Partie

 

Introduction

Je vais vous présenter l’outil de conversion de fichier IFC vers PostGIS que je suis en train de réaliser, outil développé sous en VB.Net  2008 avec :

• le composant IFCsvr ActiveX de SECOM Co. Ltd.,
• la librairie IFC Engine DLL de TNO Building Research .

Cet outil doit permettre d’exporter les données d’un fichier IFC vers une base de données PostGIS, ce billet s’enrichira au fur et à mesure de l’évolution du développement de l’outil.

IFCsvr ActiveX

Le composant ActiveX IFCsvr est un composant ActiveX pour le traitement de données IFC(Industry Foundation Classes) . Ce composant fournit un environnement de programmation compatible IFC avec Microsoft Visual Basic, Visual Basic pour Application, VC + +, C #, Visual Basic.Net, ...

Ce composant est gratuit pour un usage non commercial.

Voici le lien pour le groupe ifcsvr-users : ifcsvr-users · IFCsvr ActiveX Component

Il vous faudra vous référencer sur le groupe pour pouvoir télécharger le composant (IFCsvrR300) et autres documents.

IFC Engine DLL

Voici les lien pour la librairie IFC Engine : IFCEngine ou IFCdll

Cette librairie est gratuite pour un usage non commercial.

 

 

Industry Foundation Classes

 

Ci-dessous des extraits de l’article Wikipédia sur le format IFC :

Source : Article Industry Foundation Classes de Wikipédia en français (http://fr.wikipedia.org/wiki/Industry_Foundation_Classes).

Début citation.

“Le format IFC (Industry Foundation Classes) est un format de fichier orienté objet utilisé par l'industrie du bâtiment pour échanger et partager des informations entre logiciels.

Origine des IFC

Ce standard est né de l'initiative de l'IAI (International Alliance for Interoperability), associant des entreprises du secteur de la construction et des éditeurs de logiciels. Cette organisation est divisée en chapitres nationaux dont BuildingSmart est le représentant français. Elle a pour objectif principal de favoriser l'interopérabilité des logiciels dans le secteur de la Construction.

Building Information Modeling

IFC est le format usé communément par le BIM (Building Information Modeling)

Dans l'acceptation du BIM, le format IFC, standard ouvert compréhensible de tous les acteurs de la construction, organise les objets de l'industrie de la construction, les familles, autour d'un modèle informatique 3D, le BIM . Ces objets supportent le cycle de vie complet d'un bâtiment incluant les débuts-la conception, la documentation et la construction-, l'exploitation des bâtiments, la gestion des installations (facility management) et, enfin, la démolition et l'élimination².

Spécification des IFC

Le modèle de données IFC est spécifié en langage EXPRESS, conforme à la norme ISO 10303-11 (STEP part 11). C'est un modèle orienté objet qui définit des classes associées à tous les objets de construction. Un modèle d'ouvrage est donc un ensemble hiérarchisé d'instances de classes du modèle IFC. Les échanges s'effectuent par l'intermédiaire de fichiers STEP-21 (ISO 10303-21) ou de fichiers XML (ifcXML).

Par exemple, la hiérarchisation d'un bâtiment s'effectue ainsi :

• une zone contient plusieurs bâtiments.
• un bâtiment contient plusieurs étages
• un étage contient plusieurs salles etc.

…

Formats

Formats IFC⁴

Format	Définition
.ifc	Fichier IFC utilisant la structure STEP conformément à ISO10303-21. Structuration des données selon les spécifications Express. C'est le format d'échange par défaut
.ifcXML	Fichier IFC utilisant la structure XML. Le format .ifcXML peut être généré directement par l'application ou à partir d'un fichier .ifc, conformément à ISO10303-28, la représentation XML des schémas et données Express. Le format .ifcXML est normalement 300-400% plus gros que le fichier .ifc
.ifcZIP	Fichier IFC utilisant l'algorithm de compression PKzip 2.04g

…

Lien externe

• Site du chapitre francophone de l'IAI (International Alliance for Interoperability)
• Site technique de l'IAI
• Wiki dédié au ifc

• Site du Projet Expert En france, le projet Expert rassemble les principaux acteurs du monde de la construction et du Bâtiment, il vise à promouvoir l'utilisation des IFC et de la BIM dans ce secteur. Il est soutenu par l'État (ministère chargé de la construction).

“

Fin citation.

 

Installation

Téléchargez l’outil IFC2PG : ifcpg.rar

Décompressez l’archive sous c:\

Créez un raccourci sur le bureau.

 

Utilisation

 

Lancez l’application par un double-clic sur le raccourci.

 

L’interface s’enrichira au fur a mesure du développement, la barre d’outil principal ne comporte que deux icônes :

Cet icône permet d’ouvrir un fichier IFC.

Cet icône permet de se connecter à une base de données PostgreSQL/PostGIS.

Ouvrez le fichier qui se trouve dans le répertoire IFC

1. Arborescence du fichier IFC, les icônes avec une ampoule peuvent être éteints (dé affichage de tous les objets enfants) ou allumés (affichage de tous les objets enfants)
2. Zone d’affichage du fichier IFC (IFCEngine OCX)
3. Barre d’outils pour interagir avec l’affichage :

                

1. Rotation 3D
2. Zoom
3. Déplacement
4. Vue de Face
5. Vue de Coté

Dé afficher les deux derniers IfcBuildingStorey

Connecter vous à une base PostGIS

Export des données vers PostGIS

A ce stade du développement, pour réaliser un export il faut faire un clic droit sur un objet pour ouvrir un menu qui contiendra l’ordre d’export. Seul les objets IfcWallStandardCase (Mur) peuvent être exportés, un par un en sélectionnant chaque IfcWallStandardCase  ou bien en sélectionnant l’objet IfcBuildingStorey (Etage) pour exporter tous les objets IfcWallStandarCase.

Lorsque vous cliquez sur un objet appartenant à un objet IfcBuildingStorey, celui-ci passe en surbrillance.

Sélectionnez le 4eme objet IfcBuildingStorey, faites un clic droit et sélectionnez “Export de tous les murs” dans le menu.

La fonction va créer la table IfcWallStandardCase si elle n’existe pas, dans la base de données PostGIS à laquelle vous êtes connecté, il n’y a pas encore de gestion de doublons. La structure de la table pourra changer avec les prochaines version de l’outil. Faites moi savoir si vous avez des suggestions quand à son contenu.

 

Il ne vous reste plus qu’a allez visualiser le contenu de la table avec par exemple OpenJump.

Le champ géométrique de la table contient les axes des murs, un second champ géométrique pourrait contenir le contour, a voir.

Conclusion

Nous voila à la fin de cette première partie, l’outil réalise l’export des axes des murs (IfcWallstandardCase), je vais rajouter petit à petit l’export des autres objets.

J’attends vos remarques, suggestions ce premier jet d’un outil permettant d’exporter les objets d’un fichier IFC vers une bas de données PostGIS.

AutoHotKey, OpenJump, PgAdmin3,PostGIS et Google Earth…

Présentation

AutoHotKeyest un utilitaire open-source pour Windows. Avec lequel, vous pouvez :

• Automatiser presque n'importe quelle tâche en envoyant des frappes clavier et des clics souris. Vous pouvez écrire des macros à la main ou utiliser l’enregistreur de macro.
• Créez des raccourcis pour le clavier et la souris. Pratiquement n'importe quelle touche, bouton, ou combinaison peuvent devenir une touche directe.
• Réaliser des substitutions de chaînes de caractères, par exemple en tapant «cad» vous pouvez automatiquement obtenir « c’est à dire».
• Créer des formes de saisie de données, des interfaces utilisateurs et des barres de menu faites sur commande.
• Remappage des touches et des boutons sur votre clavier, et souris.
• Convertir n’importe quelle macro en un exécutable pouvant fonctionner sur un ordinateur ne possédant pas AutoHotKey.

Vous trouverez plus d’informations sur le site AutoHotKey au lien suivant :
http://www.autohotkey.com

Installation

Téléchargez la version 1.0.48.05 de AutoHotKey à l’adresse suivante : http://www.autohotkey.com/download/AutoHotkeyInstall.exe
Téléchargez SciTE4autoHotKey Version 3 ; l’editeur gratuit de script AutoHotKey à l’adresse suivante :
http://www.autohotkey.net/~fincs/SciTE4AutoHotkey_3/BetaInstall_v3_beta4.exe

Introduction

Voici les version des outils que nous utiliserons dans ce tutorial :

• AutoHotKey 1.0.48.05 ,
• SciTE4AutoHotKey Version 3 ,
• OpenJump 1.4.0.3,
• PgAdmin III 1.12.1 ,
• PostgreSQL 9.0.1,
• PostGIS 1.5.2

Je vais vous présenter des scripts dont vous pourrez tirer profit dans l’utilisation quotidienne de vos bases de données PostGIS.

Script 1 : Lancer OpenJump, PgAdmin III par une combinaison de touches

Ce premier script une fois compilé permettra de lancer :

• OpenJump 1.4.0.3avec la combinaison de touche suivantes :

+

• PgAdmin III avec la combinaison de touches suivantes :

+

Ouvrir SciTE4AutoHotkey (C:\Program Files\AutoHotkey\SciTE_beta4\scite.exe), copier le texte ci dessous dans l’éditeur.

#NoEnv          ; Recommended for performance and compatibility with future AutoHotkey releases.
#Persistent
#SingleInstance
SendMode Input  ; Recommended for new scripts due to its superior speed and reliability.
;--- touche Windows+o  : lance OpenJump ----
#o::
run C:\Program Files\OpenJUMP1.4.0.3\bin\bin\OpenJUMP.exe
return
;--- touche Windows+p  : lance PgAdmin III ---
#p::
run,C:\Program Files\PostgreSQL\9.0\bin\pgAdmin3.exe
return

Sauvegarder le fichier sous AH_ToolsExec.ahk, puis appuyer sur la touche F5 pour lancer l’exécution.
lien pour télécharger AH_ToolsExec.ahk


Vous devriez voir apparaitre un icône en forme de H dans la barre de tâches.

L’appui des séquences de touches définis dans le script lancera son programme associé.
Pour arrêter le script, faites un clic droit sur l’icône, puis choisissez exit dans le menu.

Nous avons vu comment lancer l’exécution de programmes à partir de combinaisons de touches, vous pouvez enrichir ce script en rajoutant d’autres programmes.

Script 2 : Piloter OpenJump

Nous allons créer le script AH_OpenJump.ahk qui va permettre de piloter OpenJump par le biais de combinaison de touches, mais c’est un des intérêts de ce script, qui ne seront  actives que si la présence d’OpenJump est détecté.

Touche Windows + Touche F1

+

Cette combinaison de touche va permettre d’ouvrir la boite “Exécuter une requête SQL” et définir une Connexion, il ne vous restera plus qu’a taper une requête puis l’exécuter.

La première chose à vérifier, est que vous disposez d’une connexion. Pour cela sous OpenJump, allez dans le menu Couche, puis Exécuter une requête SQL, puis cliquez sur l’icône Gestionnaire de connexion.

Le gestionnaire de connexion s’affiche, si la liste est vide, cliquez sur le bouton Ajouter.

Créez alors une connexion en remplissant tous les champs, puis cliquez sur le bouton OK. Le nom que vous allez choisir sera utilisé dans le script, il est peut être différent du champ Database. Vous pouvez créer autant de connexion que vous avez de bases de données.

Cliquez sur le bouton OK dans le Gestionnaire de connexion, puis fermez “Exécuter une requête SQL”

Ouvrir SciTE4AutoHotkey (C:\Program Files\AutoHotkey\SciTE_beta4\scite.exe), copier le texte ci dessous dans l’éditeur.

#NoEnv          ; Recommended for performance and compatibility with future AutoHotkey releases.
#Persistent
#SingleInstance
SendMode Input  ; Recommended for new scripts due to its superior speed and reliability.
SetTitleMatchMode,1
DetectHiddenWindows,on
SetTimer, checkOpenJump, 250
return
checkOpenJump:
OpenJumpID := WinExist("OpenJUMP")
If OpenJumpID > 0
{
   ;----------------------------------------------------------------------------
    ;--- ces combinaisons de touche ne fonctionne que si OpenJump est ouvert  ---
    ;----------------------------------------------------------------------------
   
    ;--- touche Windows+F1 : Ouvre la Fenêtre "Exécuter une requête SQL"      ---
    #F1::
        IfWinNotActive, OpenJUMP, , WinActivate, OpenJUMP,
        WinWaitActive, OpenJUMP,
        ;WinActivate, ahk_id %OpenJumpID%
        ;Sleep 100
       
        WinGetPos, X, Y, Width, Height
        ToolTip, Ouverture de la fenêtre "Exécuter une requête SQL",(%X%+%Width%)/2,(%Y%+%Heigth%)/2
        SetTimer, RemoveToolTip, 5000
       
        ;--- menu Couche ---
        send !c       
        sleep 30
       
       ;--- Exécuter une requête SQL... ---
        send {Down}
        send e     
        ;--- attente ouverture fenêtre ---
        WinWaitActive,Exécuter une requête SQL
       
       ;--- Connexion ---
        MouseClick, left, 528, 39
        sleep 30
        send {Esc}
        sleep 30
        send MYBD ;--- Le nom que vous avez saisi dans ajouter une connexion
        sleep 30
       
        ;--- Requête ---
        MouseClick, left, 201, 122
        sleep 30
        WinGetPos, X, Y, Width, Height
        ToolTip, Tapez votre requête,(%X%+%Width%)/2,(%Y%+%Heigth%)/2
        SetTimer, RemoveToolTip, 5000
    return
   ;--- insérer ci-dessous la suite du script ---
}
return
RemoveToolTip:
SetTimer, RemoveToolTip, Off
ToolTip
return

Sauvegarder le fichier sous AH_OpenJump.ahk, puis appuyer sur la touche F5 pour lancer l’exécution, puis appuyez simultanément sur la touche Windows et la touche F1, cela doit déclencher l’ouverture de la boite de dialogue “Exécuter une requête SQL”

Touche Windows + Touche F2

+

Cette combinaison de touche va permettre d’ouvrir la boite “Exécuter une requête SQL” ,définir une Connexion, puis exécuter toutes les requêtes contenues. Le fichier de requêtes (requetes_OpenJump.txt )doit se trouver dans le même répertoire que le script.
Insérer la suite du script dans le fichier AH_OpenJump.ahk.

    ;--- touche Windows+F2 : Ouvre la Fenêtre "Exécuter une requête SQL"                              ---
    ;                                     Exécute chaque requête contenue dans le fichier requetes.txt,     ---
    ;                                     chaque requête doit être séparée par ";"
    #F2::
       ;--- vide le presse papier --
        clipboard =
       
       ;--- lit le fichier requête ---
        ligne=
        tabRequetes=
        sep:=";" ;--- séparateur de requêtes
       
        Loop, Read, requetes_OpenJump.txt ;--- le fichier de requêtes
        {
            ligne = %A_LoopReadLine%
            tabRequetes = %tabRequetes%%ligne% `n
            Sleep,100

            IfInString, tabRequetes, %sep%
            {
                clipboard = %clipboard%%tabRequetes% `n
                ClipWait
                Sleep,500
                               
               ;------------------
                WinActivate, ahk_id %OpenJumpID%
                Sleep 100
               
                ;--- menu Couche ---
                send !c       
                sleep 30
               
                ;--- Exécuter une requête SQL... ---
                send {Down}
                send e     
               ;--- attente ouverture fenêtre ---
                WinWaitActive,Exécuter une requête SQL
               
                ;--- Connexion ---
                MouseClick, left, 528, 39
                sleep 30
                send {Esc}
                sleep 30
                send MYBD;--- Le nom que vous avez saisi dans ajouter une connexion
                sleep 30
               
                ;--- Requête ---
                MouseClick, left, 201, 122
                sleep 30               
               
               ;--- si l'éditeur de requêtes est présent ---
                queryID := WinExist("Exécuter une requête SQL")
                if queryID > 0
                {
                    WinActivate
                   ;--- colle le contenu du presse papier dans le requêteur ---
                    send ^a
                    sleep 100
                    send ^v
                    Sleep,100
                   
                    ;--- bouton OK ---
                    send !o
                    sleep 100
                }
                clipboard =
                tabRequetes=
                ligne=
                sleep 100
            }
        }
       
       
    return

Créez le fichier requetes_OpenJump.txt, avec une ou plusieurs requêtes, ci-dessous un exemple. Il y a deux contraintes, utiliser la fonction ST_Asbinary pour les géométries et séparer chaque requête par un point virgule.

SELECT st_asbinary(wkb_geometry) as geom,ogc_fid,nom
FROM commune;
SELECT st_asbinary(wkb_geometry) as geom,ogc_fid
FROM batiment;

Sauvegarder le fichier sous AH_OpenJump.ahk, puis appuyer sur la touche F5 pour lancer l’exécution, puis appuyez simultanément sur la touche Windows et la touche F2, cela doit déclencher l’ouverture de la boite de dialogue “Exécuter une requête SQL”, choisir une connexion, puis exécuter les requêtes que vous aurez écrites.

voici le lien pour télécharger AH_OpenJump.ahk

Script 3 : Exporter le résultat d’une requête dans un fichier KML et l’ouvrir sous Google Earth

Nous allons créer le script AH_PG2GE.ahk qui va permettre d’exécuter une requête contenue dans un fichier, exporter le résultat de la requête dans un fichier KML et l’ouvrir sous Google Earth.
Pour pouvoir utiliser ce script, vous devez avoir un pilote ODBC pour PostgreSQL d’installé, si ce n’est pas le cas, voici un lien pour le télécharger :
psqlodbc_09_00_0200

Touche Windows + Touche g

+

Ouvrir SciTE4AutoHotkey (C:\Program Files\AutoHotkey\SciTE_beta4\scite.exe), copier le texte ci dessous dans l’éditeur.
#NoEnv          ; Recommended for performance and compatibility with future AutoHotkey releases.
#Persistent
#SingleInstance
#Include Com.ahk
SendMode Input  ; Recommended for new scripts due to its superior speed and reliability.
   
;--- touche Windows+g :
;--- execute le contenu du fichier exportkml.txt , cree le fichier kml et l'ouvre dans googleearth
#g::   
    COM_init()
    sQuery =
    ;--- lit le fichier de requêtes ---
    Loop, Read, exportkml.txt
    {
        ligne = %A_LoopReadLine%
        sQuery = %sQuery%%ligne%`n
    }   
    if strlen(sQuery)>0
    {
        ;--- Chaine de connection ---
        sConnect := "Driver={PostgreSQL ANSI};Server=localhost;Port=5432;Database=MYBD;Uid=postgres;Pwd=****;"
        sConnect.=  "cache Size=100, Max Varchar=254, Max LongVarChar=-4;Extended Properties=""MaxVarcharSize=1000000;MaxLongVarcharSize=300000"""
       
        ;---- cree un objet recordset que va remplir la requête ---
        prs := COM_CreateObject("ADODB.Recordset")
        COM_Invoke(prs, "Open", sQuery, sConnect)
       
       ;--- suppression fichier ---
        FileDelete,result.kml
        ;--- ouverture fichier kml ---
        FileAppend, <?xml version="1.0" encoding="utf-8" ?>`n, result.kml
        FileAppend, <kml xmlns="http://earth.google.com/kml/2.1">`n, result.kml
        FileAppend, <Document><Folder><name>result</name>`n, result.kml
       
       ;--- lecture recordset ---
        Loop
        {
            sData =
            isData:=0
           If   COM_Invoke(prs, "EOF")
              Break
             
            FileAppend, <Placemark>`n, result.kml
            FileAppend, <description><![CDATA[`n, result.kml

           ;--- lecture des champs de l'enregistrement courant ---
            pFields := COM_Invoke(prs, "Fields")
            nNumFields := COM_Invoke(pFields, "Count")
            Loop, % nNumFields
            {
                pField := COM_Invoke(pFields, "Item", A_Index-1)
                sName := COM_Invoke(pField, "Name")
                sName = %sName%
                sData := COM_Invoke(pField, "Value")
                sData = %sData%
               
                IfInString,sData, <coordinates>
                {
                    if isData = 1
                        FileAppend, ]]></description>`n, result.kml
                    FileAppend, %sData%`n, result.kml
                }
                else
                {
                    FileAppend,<b> %sName% :</b> <i> %sData%</i><br />`n, result.kml
                    isData:=1
                }
                COM_Release(pField)
            }
            COM_Release(pFields)
           
            FileAppend, </Placemark>`n, result.kml       
           
            ;--- on passe a l'enregistrement suivant ---
            COM_Invoke(prs, "MoveNext")             
        }
        ;--- fermeture recordset ---
        COM_Invoke(prs, "Close")
        COM_Release(prs)
        ;--- fermeture fichier kml ---
        FileAppend, </Folder></Document></kml>`n, result.kml
       
        ;--- ouverture dans Google Earth ---
        run,result.kml       
    }
    COM_Term()
return

Créez le fichier exportkml.txt, avec une requête, ci-dessous un exemple :

SELECT nom,ST_Askml(2,wkb_geometry,6)
FROM commune;

Sauvegarder le fichier sous AH_PG2GE.ahk, puis appuyer sur la touche F5 pour lancer l’exécution, puis appuyez simultanément sur la touche Windows et la touche g.

Vous devriez voir apparaitre le résultat sous Google Earth, le fichier kml est créé dans le répertoire du script.

Conclusion

Voilà, j’espère que ces quelques scripts vous auront donné envie d’aller plus loin, allez sur le site d’AutoHotKey, consultez le forum, épluchez la documentation, créez vos propres scripts.
Et pourquoi ne pas lancer un concours de scripts, ici ou ailleurs, j’attend vos propositions.
a vos claviers, prêt, partez….

Squelettisation

Du diagramme de Voronoï à la Squelettisation de polygone

Présentation

Boston GIS (pour Boston Geographic Information Systems http://www.bostongis.com/) est un excellent site que je consulte régulièrement. Il y a environ deux ans, j'ai lu une série de trois tutoriels présentant PL/R. Ces lectures ont fait germer en moi l'idée de la mise en place d'un algorithme de squelettisation de cours d'eau surfacique (POLYGON) dans l'environnement PostGIS.

J'ai mis au point cet algorithme il y a un an, et très récemment en téléchargeant l'extension Skeletonizer 1.0 pour OpenJump , j'ai eu la surprise après la lecture du guide utilisateur de remarquer une approche similaire pour la réalisation de cette fonctionnalité de squelettisation.


Je vais donc vous présenter la mise en place de l'algorithme tel que je l'ai défini, mais je vous invite à étudier celui défini dans le guide utilisateur de Skeletonizer, en sachant que l'équipe de développement est issue de Refractions Research , a qui l'on doit PostGIS, uDig, Geotools, Geoserver.


PostgreSQL/PostGIS quelles versions ?

La version de PostgreSQL que j'utilise est 9.0.1 (32 bit).
La version de PostGIS que j'utilise est la 1.5.2.3 (32 bit).


Qu'est-ce que la Squelettisation ?

La squelettisation consiste à réduire une forme en un ensemble de courbes, appelées squelettes, centrées dans la forme d'origine.


Qu'est-ce qu'est PL/R ?

PL / R est une extension du langage PostgreSQL qui permet d'écrire des fonctions PostgreSQL et des fonctions d'agrégation dans le langage de calcul statistique R.


Qu'est-ce que le langage R ?

R est à la fois un langage ainsi qu'un environnement pour faire de l'analyse statistique. R est disponible en logiciel libre sous licence GPL. Pour ceux qui connaissent des environnements tels que S, MatLab, et SAS - R a le même objectif. Pour plus d'informations sur R consulter le lien suivant :

http://www.r-project.org/

Lien pour chaque tutorial de BostonGIS :

PL/R Part 1
PL/R partie 1 (traduction automatique par Google)
PL/R Part 2 : PL/R and PostGIS
PL/R partie 2 (traduction automatique par Google)
PL/R Part 3
PL/R partie 3 (traduction automatique par Google)

Je vous recommande la lecture de chaque tutorial, mais plus particulièrement :

• le premier (PL/R Part 1) qui vous permettra d'installer R, ainsi que l'extension PL/R.
  
• le second (PL/R Part 2) qui vous permettra d'installer le package deldir pour R, qui est indispensable pour la réalisation de ce qui va suivre.
  

Qu'est ce que permet deldir ?

Ce package permet de réaliser :

• Une triangulation de Delaunay,
  
• Un diagramme de Voronoï. C'est cette fonctionnalité que nous allons utiliser.
  

Qu'est-ce qu'un diagramme de Voronoï ?

Un diagramme de voronoï forme une partition de l'espace qui attribue à chaque point une parcelle englobante qui comprend tout l'espace où l'on se trouve plus proche de ce point que de tout autre. Il doit son nom au mathématicien russe Georgi Fedoseevich Voronoï (1868 - 1908).

La figure ci-dessous montre la construction d'un diagramme de Voronoï à partir d'un ensemble de points.

Introduction

Voila le décor est planté nous allons pouvoir passer aux choses sérieuses, dans un premier temps nous allons créer un schéma sous PostgreSQL qui contiendra les fonctions et les tables nécessaire à la réalisation de ce tutoriel. Nous utiliserons PgAdmin III 1.12.1 qui est compatible avec la version de PostgreSQL 9.0, pour la connexion à la base de données et la réalisation de toutes les requêtes.

A la fin de chaque étape de création d'objets géométriques nous visualiserons le résultat sous OpenJump (à télécharger impérativement) dont la dernière version est compatible avec PostgreSQL 9.0. A la fin de ce tutoriel nous verrons comment utiliser le plugin PostGisViewer que j'ai développé pour PgAdmin III, pour réaliser l'enchainement de toutes les requêtes et en visualiser le résultat. Plugin qui a été présenté à ma grande satisfaction dans le Posgres on line journal en même temps que la version 1.13 de PgAdmin III (encore en phase de développement) qui apporte entre autre une nouvelle gestion des plugins, vous pouvez consultez cet article à l'adresse suivante :

http://www.postgresonline.com/journal/archives/180-pgAdmin113plugins_postgis.html


Etape 1 : Création du Schéma

Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III et exécutez la

-- Schema: skeleton
-- DROP SCHEMA skeleton;

CREATE SCHEMA skeleton

AUTHORIZATION postgres;
GRANT ALL ON SCHEMA skeleton TO postgres;
GRANT ALL ON SCHEMA skeleton TO public;



Etape 2 : Création de la table avec insertion d'une géométrie

Téléchargez la requête suivante et exécutez la dans l'éditeur de requêtes de PgAdmin III
 CREATE TABLE skeleton.water_area...
 
Visualisation du contenu de la table

 Nous allons visualiser le contenu de notre table sous OpenJump, pour cela allez dans le menu Couche

  Puis exécuter la requête suivante :

Et voilà à quoi ressemble le polygone que nous allons squelettiser, c'est une partie de cours d'eau avec deux iles, donc nous avons un POLYGON composé de deux « interioring » et d'un « exteriorring ».

Etape 3 : Création des fonctions


Nous allons créer maintenant les fonctions nécessaires au traitement de squelettisation.

1)  skeleton.vronoi

Cette première fonction est adaptée de la fonction voronoi  tirée du tutorial PL/R Part 2 : PL/R and PostGIS  de BostonGis.

Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III et exécutez la

/*--- fonction pour créer le diagramme de voronoi ---*/
-- Function: skeleton.voronoi(text, text, text)
CREATE OR REPLACE FUNCTION skeleton.voronoi(text, text, text)
RETURNS void AS
$BODY$
library(deldir)
# select the point x/y coordinates into a data frame...
points <- pg.spi.exec(sprintf("select x(%2$s) as x, y(%2$s) as y from %1$s;",arg1,arg2))

voro = deldir(points$x, points$y, digits=22, frac=0.00000000000000000000000001,list(ndx=2,ndy=2,fctr=0.000001), rw=c(min(points$x)-abs(min(points$x)-max(points$x)), max(points$x)+abs(min(points$x)-max(points$x)), min(points$y)-abs(min(points$y)-max(points$y)), max(points$y)+abs(min(points$y)-max(points$y))))

# dirsgs
for (i in 1:length(voro$dirsgs[[1]]))
{
#--- voronoi ---
pg.spi.exec(sprintf("INSERT INTO skeleton.%1$s (id,idedge,wkb_geometry) VALUES (%2$d,%2$d,GeomFromtext('LINESTRING(%3$f %4$f, %5$f %6$f)',4326));",
arg3,i,
voro$dirsgs[i,1],voro$dirsgs[i,2],
voro$dirsgs[i,3],voro$dirsgs[i,4]))
}
$BODY$
LANGUAGE plr VOLATILE
COST 100;
ALTER FUNCTION skeleton.voronoi(text, text, text) OWNER TO postgres;
GRANT EXECUTE ON FUNCTION skeleton.voronoi(text, text, text) TO postgres;
GRANT EXECUTE ON FUNCTION skeleton.voronoi(text, text, text) TO public;

2) skeleton.sridUtmByLongitude

  Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III et exécutez la

/*--- fonction pour obtenir le srid UTM zone de la geometrie ---*/
-- Function: skeleton.sridUtmByLongitude(geometry)
CREATE OR REPLACE FUNCTION skeleton.sridUtmByLongitude(geometry)
RETURNS integer AS
$BODY$
DECLARE
gg ALIAS FOR $1; -- Geometrie
pt geometry; -- Geometrie de type POINT
sr int4; -- srid

BEGIN
sr := 0;

IF GeometryType(gg)='POINT' THEN
pt := gg;
ELSE
pt := ST_Centroid(ST_Envelope(gg));
END IF;

sr := SRID(pt);

IF sr = 4326 THEN
SELECT srid into sr FROM spatial_ref_sys
WHERE srtext like
(
SELECT
CASE
WHEN int2(y(pt)) >= 0
THEN '%WGS 84 / UTM zone '||int2((x(pt)+180)/6+1)||'N%'::text
ELSE '%WGS 84 / UTM zone '||int2((x(pt)+180)/6+1)||'S%'
END
AS utmzone
);
END IF;

RETURN sr;
END;
$BODY$
LANGUAGE plpgsql VOLATILE STRICT
COST 100;
ALTER FUNCTION skeleton.sridUtmByLongitude(geometry) OWNER TO postgres;
GRANT EXECUTE ON FUNCTION skeleton.sridUtmByLongitude(geometry) TO postgres;
GRANT EXECUTE ON FUNCTION skeleton.sridUtmByLongitude(geometry) TO public;


3) skeleton.smoothvoronoiline

 Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III et exécutez la

/*--- fonction lissage ---*/
-- Function: skeleton.smoothvoronoiline(geometry)
-- DROP FUNCTION skeleton.smoothvoronoiline(geometry);
CREATE OR REPLACE FUNCTION skeleton.smoothvoronoiline(geometry)
RETURNS geometry AS
$BODY$
DECLARE
    geom ALIAS FOR $1; -- geometry
    sql1 text;
    rec1 RECORD;

BEGIN
    IF ST_NumPoints(geom)>2 THEN
        sql1:='SELECT ST_LineFromMultiPoint(ST_Collect(a.geom )) as geom
        FROM
        (
        SELECT ST_PointN(geometry('||quote_literal(geom::text)||'),generate_series(1,ST_NumPoints(geometry('||quote_literal(geom::text)||')),2)) as geom
        ) a;';
        FOR rec1 IN EXECUTE sql1 LOOP END LOOP;
        IF FOUND THEN
            RETURN ST_AddPoint(rec1.geom,ST_EndPoint(geom));
        END IF;
    ELSE
        RETURN geom;
    END IF;

END;
$BODY$
LANGUAGE plpgsql VOLATILE STRICT
COST 100;
ALTER FUNCTION skeleton.smoothvoronoiline(geometry) OWNER TO postgres;
GRANT EXECUTE ON FUNCTION skeleton.smoothvoronoiline(geometry) TO postgres;
GRANT EXECUTE ON FUNCTION skeleton.smoothvoronoiline(geometry) TO public;

Etape 4 : Squelettisation

1)  Création du diagramme de voronoï

Nous allons utiliser la requête suivante pour générer un diagramme de voronoï, pour cela nous utilisons la fonction skeleton.voronoi créée précédemment qui reçoit deux paramètres :

Paramètre 1 : une liste de points que l'on extrait de l'exteriorring et de chaque interiorring,

SELECT ST_PointN(wkb_geometry,generate_series(1, ST_NumPoints(wkb_geometry))) as geom

FROM

(

SELECT st_exteriorring(wkb_geometry) as wkb_geometry,

ST_Length(ST_Transform(st_exteriorring(wkb_geometry),skeleton.sridUtmByLongitude(wkb_geometry))) as length

FROM skeleton.water_area

) as a

UNION ALL

SELECT ST_PointN(wkb_geometry,generate_series(1, ST_NumPoints(wkb_geometry))) as geom

FROM

(

SELECT ST_InteriorRingN(wkb_geometry,generate_series(1,ST_NumInteriorRings(wkb_geometry))) as wkb_geometry,

ST_Length(ST_Transform(ST_InteriorRingN(wkb_geometry,generate_series(1,ST_NumInteriorRings(wkb_geometry))),skeleton.sridUtmByLongitude(wkb_geometry))) as length

FROM skeleton.water_area WHERE id=1

) as b

Paramètre 2 : la table qui contiendra les lignes composant le diagramme de voronoï.

'b.geom','voronoi_line'

Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III et exécutez la

SELECT skeleton.voronoi('(
SELECT ST_PointN(wkb_geometry,generate_series(1, ST_NumPoints(wkb_geometry))) as geom
FROM
(
SELECT st_exteriorring(wkb_geometry) as wkb_geometry,
ST_Length(ST_Transform(st_exteriorring(wkb_geometry),skeleton.sridUtmByLongitude(wkb_geometry))) as length
FROM skeleton.water_area
) as a
UNION ALL
SELECT ST_PointN(wkb_geometry,generate_series(1, ST_NumPoints(wkb_geometry))) as geom
FROM
(
SELECT ST_InteriorRingN(wkb_geometry,generate_series(1,ST_NumInteriorRings(wkb_geometry))) as wkb_geometry,
ST_Length(ST_Transform(ST_InteriorRingN(wkb_geometry,generate_series(1,ST_NumInteriorRings(wkb_geometry))),skeleton.sridUtmByLongitude(wkb_geometry))) as length
FROM skeleton.water_area WHERE id=1
) as b
) as b','b.geom','voronoi_line');

Visualisez le contenu la table skeleton.voronoi_line sous OpenJump, pour cela allez dans le menu Couche

Puis exécutez la requête suivante :

Et voici le resultat :

 

  Vous pouvez remarquer qu'au centre du cours d'eau se dessine un squelette, pour le rendre plus visible, exécutez la requête suivante qui va le renforcer.

Le squelette qui apparait est très irrégulier, en forme de dent de scie, et disparait dans certaines zones, plus précisément autour des deux iles. Si l'on zoom sur l'ile la plus haute, voici ce que l'on peut voir :

Le squelette sur le bras le plus étroit n'est pas présent, alors que sur le bras plus large celui est présent. A ce stade du tutoriel, il semblerait que la squelettisation de ce cours d'eau n'est pas sur la bonne voie….

Pourtant rien n'est perdu, en analysant la figure ci-dessus, il apparait que la densité de points n'est pas assez importante, PostGIS va nous fournir la fonction qui va permettre d'augmenter cette densité, et plus précisément de placer des points sur le contour extérieur ainsi que les contours intérieurs avec un espacement déterminé.

Cette fonction est la fonction ST_Line_Interpolate dont voici la définition (extraite du site www.postgis.fr ) :

Line_Interpolate(LINESTRING, location)

Retourne un point interpolé le long d'une ligne. Le premier argument doit être une LINESTRING. Le second argument est de type float8 dont la valeur doit être comprise entre 0 et 1 et qui représente une fraction de la longueur 2d totale où le point doit être localisé.

Nous allons placer des points tout les 10m le long de chaque contour.


Copiez/Collez les requêtes suivante dans l'éditeur de requêtes de PgAdmin III et exécutez les
TRUNCATE TABLE skeleton.voronoi_line;

SELECT skeleton.voronoi('(
SELECT ST_Transform(ST_Line_Interpolate_Point(wkb_geometry,generate_series(10, floor(length)::int, 10)/length),4326) as geom
FROM
(
SELECT ST_Transform(st_exteriorring(wkb_geometry),skeleton.sridUtmByLongitude(wkb_geometry)) as wkb_geometry,
ST_Length(ST_Transform(st_exteriorring(wkb_geometry),skeleton.sridUtmByLongitude(wkb_geometry))) as length
FROM skeleton.water_area
) as a
UNION ALL
SELECT ST_Transform(ST_Line_Interpolate_Point(wkb_geometry,generate_series(10, floor(length)::int, 10)/length),4326) as geom
FROM
(
SELECT ST_Transform(ST_InteriorRingN(wkb_geometry,generate_series(1,ST_NumInteriorRings(wkb_geometry))),skeleton.sridUtmByLongitude(wkb_geometry)) as wkb_geometry,
ST_Length(ST_Transform(ST_InteriorRingN(wkb_geometry,generate_series(1,ST_NumInteriorRings(wkb_geometry))),skeleton.sridUtmByLongitude(wkb_geometry))) as length
FROM skeleton.water_area WHERE id=1
) as b
) as b','b.geom','voronoi_line');

Puis exécutez la requête suivante pour en visualiser le résultat :

  

  Puis enchainons la requête qui va renforcer le squelette.

Cela devient intéressant, le squelette apparait très nettement que ce soit dans le bras le plus étroit comme dans la partie la plus large. Nous allons zoomer sur le bras le plus étroit pour mieux se rendre compte du résultat, en désactivant auparavant l'affichage du diagramme de voronoï.

Nous pouvons constater, que le squelette apparait légèrement en dent de scie, l'application d'une fonction de lissage semble nécessaire, nous allons donc utiliser la fonction skeleton.smoothvoronoiline pour réaliser ce lissage.


Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III et exécutez la


DROP TABLE IF EXISTS skeleton.skeleton_line;


CREATE TABLE skeleton.skeleton_line WITH OIDS AS
SELECT a.id,a.geom,st_Length_Spheroid(a.geom,'SPHEROID("WGS 84",6378137,298.257223563)') as longueur,ST_NumPoints(a.geom) as nbpt
FROM
(
SELECT (st_dump(ST_LineMerge(ST_Union(wkb_geometry)))).path[1] as id,skeleton.smoothvoronoiline((st_dump(ST_LineMerge(ST_Union(wkb_geometry)))).geom) as geom
FROM skeleton.voronoi_line
) a;

Puis visualisons le résultat sous OpenJump :

 

Après le passage de la fonction de lissage, nous obtenons un squelette ou l'effet dent de scie a disparu. Pour avoir une vue d'ensemble du squelette, zoomons sur l'étendue maximale du cours d'eau, regardons ci-dessous le résultat.

Sur cette vue d'ensemble, nous pouvons remarquer la présence de petites branches sur le squelette principale. Nous allons essayer de simplifier encore ce résultat, en éliminant certaines de ces branches. Affichons sous OpenJump les attributs de la table skeleton .skeleton_line qui sont triés sur le critère de la longueur de chaque LINESTRING.

 Nous remarquons tout de suite que la longueur de certaines branches est inférieure à 10m (9 branches), inférieure à 70m (45 branches). Ce critère de longueur va nous permettre d'éliminer celles qui seront en dessous de ce critère et dont une extrémité ne sera pas connectée a une autre partie du squelette.

Avant de les supprimer, nous allons afficher les branches candidates à la suppression par le biais de la requête ci-dessous, en fixant le critère à 70m.


SELECT ST_AsBinary(geom),longueur
FROM skeleton.skeleton_line
WHERE id NOT IN(
SELECT a.id
FROM
(
SELECT a.geom,a.id,a.longueur,a.nbpt,(select count(b.*)
FROM skeleton.skeleton_line b
where (Distance(startpoint(a.geom),startpoint(b.geom))=0 or Distance(startpoint(a.geom),endpoint(b.geom))=0))-1 as nbconS,
(SELECT COUNT(b.*) FROM skeleton.skeleton_line b where (Distance(endpoint(a.geom),startpoint(b.geom))=0
or distance(endpoint(a.geom),endpoint(b.geom))=0))-1 as nbconE
FROM skeleton.skeleton_line a
) a
WHERE NOT (a.longueur < 70.0 and a.nbconS=0 or a.nbconE=0)
ORDER BY a.id
);

 

Les branches qui répondent au critère de longueur et dont une extrémité n'est pas connectée, apparaissent en rouge sur la figure ci-dessus. Elles peuvent être supprimées sans que la structure du squelette ne soit remise en cause, ce que nous allons réaliser avec la requête ci-dessous.


Copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III


DELETE FROM skeleton.skeleton_line
WHERE id NOT IN(
SELECT a.id
FROM
(
SELECT a.geom,a.id,a.longueur,a.nbpt,(SELECT COUNT(b.*) FROM skeleton.skeleton_line b WHERE (Distance(startpoint(a.geom),startpoint(b.geom))=0 or Distance(startpoint(a.geom),endpoint(b.geom))=0))-1 as nbconS,
(SELECT COUNT(b.*) FROM skeleton.skeleton_line b WHERE (Distance(endpoint(a.geom),startpoint(b.geom))=0 or distance(endpoint(a.geom),endpoint(b.geom))=0))-1 as nbconE
FROM skeleton.skeleton_line a
) a
WHERE NOT (a.longueur < 70.0 AND a.nbconS=0 or a.nbconE=0)
AND a.id
);


Exécutez la jusqu'à ce qu'il n'y ait plus de ligne modifiée.

La requête a été exécutée avec succés : 32 lignes modifiées. La requête a été exécutée en 177 ms.
La requête a été exécutée avec succés : 3 lignes modifiées. La requête a été exécutée en 22 ms.
La requête a été exécutée avec succés : 0 ligne modifiée. La requête a été exécutée en 27 ms.

Puis visualisons le résultat sous OpenJump :

 
  
Voila nous y sommes, le squelette est créé, il reste quelques branches qui pourraient être supprimées en augmentant le critère de longueur, pour notre exemple nous nous arrêterons là.


    
Utilisation du plugin PostGisViewer pour PgAdmin III


Pour la mise en place du plugin je vous renvoie à ce lien :

Plugin PostGisViewer pour PgAdmin III


Une fois la mise en place du plugin effective, si vous utilisez une version de PostgreSQL dont la version est supérieur ou égale à la 9.0, copiez/Collez la requête suivante dans l'éditeur de requêtes de PgAdmin III, exécutez la, puis fermez PgAdmin IIII .

ALTER DATABASE "NOM DE VOTRE BASE" SET bytea_output='escape';

Cette astuce est extraite de l'article  « PgAdmin III 1.13 - change in plugin architecture and PostGIS Plugins »  dans la partie « PostGIS plugin for displaying geometries in PgAdmin III »  dont je vous conseille la lecture.

1. Vous pouvez maintenant rouvrir PgAdmin III, copiez/collez la requête suivante dans l'éditeur de requête (un seul éditeur ouvert à la fois)
   

SELECT st_asbinary(wkb_geometry) as geom,geometrytype(wkb_geometry),id

FROM skeleton.water_area;

SELECT ST_AsBinary(wkb_geometry) as geom,geometrytype(wkb_geometry)

FROM skeleton.voronoi_line;

SELECT ST_AsBinary(geom) as geom, geometrytype(geom),longueur

FROM skeleton.skeleton_line

ORDER BY longueur;

1. Dans le navigateur d'objets de PgAdmin, vous devez être positionné sur votre base de données.
   
2. Dans le menu plugins sélectionnez PostGisViewer, vous devriez obtenir ceci :
   

Sous Map Layers, chaque layer est le résultat d'une requête.


Conclusion


En premier lieu je tiens à remercier :

Georgi Fedoseevich Voronoï, sans qui rien de tout cela ne serait arrivé,
BostonGIS pour ses tutoriels,
Joe Conway pour PL/R, qui a rendu possible la mise en place de cet algorithme,
Postgres OnLine journal pour ses articles,
Refractions Research pour ses développements (PostGIS tout particulièrement),
...

Ce tutoriel est terminé, nous avons vu comment grâce au diagramme de voronoï, réaliser la squelettisation d'un cours d'eau de type POLYGON. L'enchainement de toutes ces requêtes pourrait bien sur être facilement remplacé par une fonction, à qui l'on passerait deux paramètres :

• la géométrie à traiter,
  
• la tolérance pour la suppression des branches.
  

A vous de jouer…