Le blog de Jean David TECHER, un Réunionnais à Saint-Priest/Lyon

Objectif

Bon ce soir je me suis lancé le défi de pouvoir enfin afficher la portion de la RN113 qui traverse le Gard (département n° 30) dans GoogleEarth.

1 - Ce qui ne marche pas!

Je dispose pour celà  de la table des tronçons! Donc pour connaà®tre la route en question il me suffit de faire

dde30=#  SELECT id_1,id_2,route,geometrytype(the_geom),numgeometries(the_geom) FROM route_dep30 WHERE route='n113' LIMIT 10;
   id_1    | id_2 | route |  geometrytype   | numgeometries
-----------+------+-------+-----------------+---------------
 300005729 |  204 | n113  | MULTILINESTRING |             1
 300005730 |  205 | n113  | MULTILINESTRING |             1
 300005775 |  208 | n113  | MULTILINESTRING |             1
 300005786 |  210 | n113  | MULTILINESTRING |             1
 300005787 |  211 | n113  | MULTILINESTRING |             1
 300005816 |  212 | n113  | MULTILINESTRING |             1
 300005823 |  213 | n113  | MULTILINESTRING |             1
 300005840 |  215 | n113  | MULTILINESTRING |             1
 300005853 |  216 | n113  | MULTILINESTRING |             1
 300005876 |  217 | n113  | MULTILINESTRING |             1
(10 lignes)

Donc fier de moi comme Djay (la semaine dernière) m'avait appris une fonction de PostGIS que je ne connaissais pas linemerge qui simplifie par exemple un MULTLINESTRING contenant uniquement un simple LINESTRING en un simple LINESTRING, je me dis ben c'est gagné!

Ben non ......!

En effet,cette table contient des MULTILINESTRINGS dont le nombre de LINESTRINGS est bien > 1 comme le prouve cette requête

dde30=#  SELECT distinct route,geometrytype(the_geom),numgeometries(the_geom) FROM route_dep30 WHERE route='n113';
 route |  geometrytype   | numgeometries
-------+-----------------+---------------
 n113  | MULTILINESTRING |             1
 n113  | MULTILINESTRING |             2
 n113  | MULTILINESTRING |             4
(3 lignes)

Le nouveau switch -S de shp2pgsql qui permet de convertir un MULTI en simple ojet spatial le confirme aussi:

postgres@jenna-edgy:~$ shp2pgsql -DIdS route_dep30.shp route_dep30|iconv -f LATIN1 -t UTF-8|psql dde30
Shapefile type: Arc
Postgis type: LINESTRING[2]
We have a MultiLineString with 2 parts, can't use -S switch! <-------- ICI !!!
... ...

Au premier abord j'avais cru comprendre que linemerge() transformait un MULTILINESTRING en un LINESTRING. En fait non, erreur de ma part, suite à  une explication avec Djay. Cette fonction fait bien plus que celà . Voyons donc quelques exemples:

La fonction linestring()

Exemple 1 avec deux lLINESTRING contigà¼es mise en relation par le POINT(5 5)

testgis=# select astext(linemerge(geometryfromtext('MULTILINESTRING((0 0,5 5),(5 5,6 10))',-1)));
          astext
--------------------------
 LINESTRING(0 0,5 5,6 10)
(1 ligne)

No soucy!

Exemple 2 deux LINESTRING discontigà¼es

testgis=# select astext(linemerge(geometryfromtext('MULTILINESTRING((0 0,5 5),(5 4,6 10))',-1)));
                astext
---------------------------------------
 MULTILINESTRING((0 0,5 5),(5 4,6 10))
(1 ligne)

Donc aucun souci! Après mettre entretenu avec Djay, j'ai finalement compris comme il me le suggérait que si on voulait travailler unqiuement avec la fonction linermerge, il faudrait faire une analyse des données: trouver les MULTILINESTRING qui ne sotn aps convertis en simple LINESTRING, les ouvrir dans le MULTILINESTRING correspondant, scinder etc...Mais mon but ici est de parvenir à  faire ça en une seule requête! Donc j'ai laissé tombé pour linemerge! Apparement ce sont avec mes données que je dosi faire!

Première conclusion: linemerge() ----------> tu sors ............!

Deuxième conclusion: shp2pgsql ----------> tu sors ............!

2 - Ce qui marche!

Après une heure de recherche intensive sur la documentation de PostGIS, et sur la mailing-list, j'ai trouvé la solution ze dump() function!

En pré-requis comme gdal 1.3.2 avec ogr2ogr ne permet pas de faire de sortie en fichier kml, il faut commencer par télécharger les sources de gdal 1.4.0 beta 1 à  http://www.gdal.org/dl/ puis

wget http://www.gdal.org/dl/gdal140beta1.zip
unzip gdal140beta1.zip
cd gdal-1.4.0beta1
 CC=gcc-3.4 CXX=g++-3.4 ./configure --prefix=/opt/gdal-1.4.0 --with-geos --with-static-proj4=/usr/local/--without-python
make
make install

NB: Ici j'ai tout installé dans le répertoire /opt/gdal-1.4.0

Puis on fait pour générer le fichier RN113.kml

/opt/gdal-1.4.0/bin/ogr2ogr -f KML /root/RN113.kml PG:'host=192.168.0.1 dbname=dde30 user=postgres' -sql \
"SELECT (dump).geom FROM (SELECTt dump(transform(the_geom,4326)) FROM route_dep30 WHERE route='n113') AS foo"

Et on matte tout ça dans Google Earth

PostGIS ---> KML: Portion de la RN113 (réunifié à  partir de ses tronçons) traversant le Gard depuis Google Earth