"Sans imagination il ne pourrait y avoir création."
Albert Jacquard.

"L'imagination est plus importante que le savoir."
Albert Einstein.

jeudi 10 mars 2011

AutoHotKey, OpenJump, PgAdmin3,PostGIS et Google Earth…

Présentation


clip_image001[7]
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 :
image
+
image
  • PgAdmin III avec la combinaison de touches suivantes :
image
+
image
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
script1_b

Vous devriez voir apparaitre un icône en forme de H dans la barre de tâches.
script1_c
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.
script1_d
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

image
+
image
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.
ExecQuerySql_2
Le gestionnaire de connexion s’affiche, si la liste est vide, cliquez sur le bouton Ajouter.
ExecQuerySql_3
ExecQuerySql_4
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.
ExecQuerySql_5
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

image
+
image
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.
result_ F2

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

image
+
image
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….Sourire

mercredi 8 décembre 2010

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…


 

mardi 6 juillet 2010

Plugin PgAdmin III : Export PostGIS




Un plugin pour PgAdmin III : Export PostGIS


Présentation 
Ce plugin d'export pour PostGIS a été développé en VB.Net 2008.
Il utilise l'utilitaire ogr2ogr qui fait partie du package d'outils Open Source FWTOOLS maintenus et mis à jour par Franck Warmerdam (d'ou les initiales FW).
Pour plus d'informations, consulter ce lien :
FWTOOLS 

Et celui-ci pour ogr2ogr :
ogr2ogr


Installation 


La version de PgAdmin III utilisée pour ce tutorial est la 1.10.3, lien ci-dessous pour la télécharger :
PgAdmin3-1.10.3


La version de FWTOOLS utilisée pour ce tutoriel est la 2.4.7, lien ci-dessous pour la télécharger :
FWTools 2.4.7



Lancer l'installation et définir le répertoire d'installation :
C:\FWTools à la place C:\Program Files\FWTools2.4.7
































Le programme d'export pour PostGIS est disponible ici :


Voir la fin du tutoriel pour la nouvelle version.


Configurer PgAdmin III

PgAdmin III 1.10 donne la possibilité d'ajouter très facilement des applications externes dans le menu plugins.
Pour activer la disponibilité de ces applications:  

  • Ouvrir PgAdmin III, dans le menu Fichier sélectionner Préférences
 












  •  Dans l'onglet Général donner le chemin de l'application PgAdmin III (pgadmin3.exe)

  
















  • Sortir de PgAdmin III,
  • Aller dans le répertoire C:\Program Files\pgAdmin III\1.10, faire une copie du fichier plugins.ini,
  • Ouvrir le fichier plugins.ini et copier les lignes suivantes a la fin du fichier :
    ;
    ;
    PostGisExport (Windows)
    ;
    [Separator]
    Title=Export PostGIS
    Command="$$PGBINDIR\PostGisExport\
    PostGisExport.exe"  "host=$$HOSTNAME" "port=$$PORT" "username=$$USERNAME" "password=$$PASSWORD" "database=$$DATABASE" "schema=$$SCHEMA"  "table=$$TABLE"
    Description=Export PostGIS
    KeyFile=$$PGBINDIR\
    PostGisExport\PostGisExport.exe
    Platform=windows
    ServerType=postgresql
    Database=Yes
    ;AppliesTo=database
    SetPassword=Yes
    •  Les variables  utilisées, contiendront respectivement :
    • $$PGBINDIR       :  répertoire dans lequel se trouve pgadmin3.exe,
    • $$HOSTNAME    :  adresse du serveur,
    • $$PORT                :  port,   
    • $$USERNAME     :  utilisateur ,
    • $$PASSWORD     :  mot de passe,
    • $$DATABASE      :  base de données,
    • $$SCHEMA          :  schéma courant,
    • $$TABLE              :  table sélectionné.

  • Enregistrer le fichier,
  • Ouvrir PgAdmin III et dans le menu plugins vous devriez voir ceci :






















Comment utiliser le plugin

Passons maintenant à l'utilisation du plugin, pour cela ouvrir PgAdmin III :
  • se connecter à un serveur, 
  • choisir une base de données, 
  • sélectionner un schéma puis une table,
  • lancer l'exécution du plugin






  • vous devriez obtenir cette fenêtre :














  • Sélectionner un format d'export parmi ceux disponibles (ESRI Shapefile est le format par défaut),
  • déployer la sortie console,



















  • Cliquer sur le bouton export pour obtenir le format d'export sélectionné.



















Un répertoire correspondant au format sélectionné est créé s'il n'existe pas, dans lequel vous trouverez le ou les fichiers créés. 
Si vous le souhaitez, vous pouvez récupérer la ligne de commande dans la sortie console.




Nouvelle Version : 1ere Mise a jour
  • J'ai ajoute le choix du système de projection en sortie,
  • La détection du chemin de FWTools, ce qui permet de gérer le cas ou l'installation n'a pas été faite sous C:\FWTools
  • Ajout du format GPX



Nouvelle Version : 2eme Mise a jour

  • Ajout du format DXF,
  • Export d'une requête, ce qui peut permettre d'exporter une vue.



Important :

 La requête doit contenir la fonction srid afin de pouvoir fixer le srid source dans la boite d'export, et se terminera par un point virgule. Ne mettre qu'une seule requête dans la boîte SQL.

Exemple :

SELECT nom,code_insee,wkb_geometry,srid(wkb_geometry) FROM commune;

Dans le cas d'une requête, dans la partie destination le nom de sortie est fixé a query1, vous pouvez bien entendu le modifier en conservant l'extension, avant de cliquer sur le bouton export.

Rappel :

Le résultat de l'export sera visible sous le répertoire .\ PostGisExport\Data\"extension"

Exemple :

\ PostGisExport\Data\SHP






Conclusion 



Si d'autres formats d'export vous intéressent (parmi ceux acceptés par ogr2ogr), faite m'en part, et je les rajouterais.