Répondre au commentaire

OGR: que la force soit avec les formats virtuels

lun 09-09-2013 Martin Laloux


	Niveau	Intermédiaire
	Logiciels utilisés	OGR
	Plateforme	Windows \| Mac \| Linux \| FreeBSD

Beaucoup de Sigistes connaissent ou du moins ont entendu parler du format virtuel de GDAL (GDAL Virtual Format Tutorial ) pour les fichiers de type raster et de ses multiples possibilités (The power of GDAL virtual formats ou Image Mosaicking with GDAL), mais il me semble que peu d'entre eux connaissent ce même format pour les couches vectorielles.

Principes:

Les caractéristiques sont exposées dans OGR Virtual Format. Tout comme celui pour les rasters, il s'agit d'un simple fichier au format XML. Ce fichier VRT (.vrt) se traite comme n'importe quel fichier vectoriel (fichier shapefile, etc.):

<OGRVRTDataSource>
    <OGRVRTLayer name="points">
        <SrcDataSource>/chemin/points.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:31370</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="x" y="y"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

Dans le cas présent, il s'agit d'un simple fichier CSV qui pourrait être ouvert avec un logiciel SIG comme QGIS plus simplement. Mais, il est possible d'aller plus loin comme:

la création directe d'un fichier 3D, non reconnu come tel par QGIS, mais bien par des logiciels comme GRASS GIS (toujours à partir d'un simple fichier CSV):

<OGRVRTDataSource>
    <OGRVRTLayer name="points">
        <SrcDataSource>/chemin/points3D.csv</SrcDataSource>
        <GeometryType>wkbPoint25D</GeometryType>
        <LayerSRS>EPSG:31370</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="x" y="y" z="z"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

le changement de système de projection de la couche (depuis la version 1.10 de GDAL):

<OGRVRTDataSource>
    <OGRVRTWarpedLayer>
        <OGRVRTLayer name="test">
            <SrcDataSource>/chemin/test.shp</SrcDataSource>
        </OGRVRTLayer>
        <TargetSRS>EPSG:4326</TargetSRS>
    </OGRVRTWarpedLayer>
</OGRVRTDataSource>

la redéfinition des noms et des formats des colonnes (depuis la version 1.7 de GDAL);
le découpage d'une couche en fonction d'un polygone (depuis la version 1.7 de GDAL, voir plus loin);
la possibilité d'utiliser des indexes spatiaux;
la possibilité de charger plusieurs couches à la fois (comme avec GDAL);
la possibilité de découper une couche en fonction des attributs;
....

De plus, tous les types de support peuvent être traités: fichiers textes, fichiers shapefiles, bases de données, base de données spatiale, avec ODBC, etc.

avec une base SQLite non spatiale:

<OGRVRTDataSource>
    <OGRVRTLayer name="Formations">
        <SrcDataSource>/chemin/test.sqlite</SrcDataSource> 
        <SrcLayer>Form</SrcLayer> 
        <GeometryType>wkbPoint</GeometryType>
            <LayerSRS>EPSG:31370</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="lng" y="lat"/> 
    </OGRVRTLayer>
</OGRVRTDataSource>

avec un fichier dbf:

<OGRVRTDataSource>
    <OGRVRTLayer name="Formations">
        <SrcDataSource>/chemin/test.dbf</SrcDataSource> 
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:31370</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="lng" y="lat"/> 
    </OGRVRTLayer>
</OGRVRTDataSource>

avec un fichier comprenant une géométrie au format WKT (voir GDAL and OGR Virtual Format, reading *.csv file and creating shapefile of lines):

<OGRVRTDataSource>
       <OGRVRTLayer name="Linestest">
           <SrcDataSource SEPARATOR="SEMICOLON">/chemin/wkt.csv</SrcDataSource>
           <SrcLayer>wkt</SrcLayer>
           <GeometryType>wkbLineString</GeometryType>
           <GeometryField encoding="WKT" field="lignes"/>
       </OGRVRTLayer>
</OGRVRTDataSource>

avec des données provenant de ArcSDE (voir Connecting to ArcSDE databases)

<OGRVRTDataSource>
       <OGRVRTLayer name="layer name in qgis"> 
           <SrcDataSource>SDE:servername,5151,,user,pwd,SCHEMA.TABLENAME,SDE.DEFAULT</SrcDataSource> 
           <SrcLayer>SCHEMA.TABLENAME</SrcLayer> 
           <GeometryType>wkbMultiPolygon</GeometryType> 
           <LayerSRS>+proj=tmerc +lat_0=39.66666666666666 +lon_0=-8.131906111111112 +k=1 +x_0=200180.598 +y_0=299913.01 +ellps=intl +units=m +towgs84=-223.237,110.193,36.649 +no_defs</LayerSRS>
       </OGRVRTLayer> 
</OGRVRTDataSource>

Je ne continuerai pas avec les possibilités basiques du format pour me focaliser sur la partie la plus intéressante, les traitements avec des requêtes SQL.

Requêtes SQL

Depuis la version 1.10 de GDAL, les commandes sont de 2 types:

OGR SQL
SQLite SQL dialect (introduit dans GDAL 1.10)

OGR SQL

En introduction, commençons par un cas simple qui pourrait être réglé par le logiciel SIG mais que nous traiterons ici avec un fichier virtuel. Il s'agit d'une carte géologique dont je veux extraire le tracé d'une formation (BUR ici) particulière. J'utiliserai QGIS pour visualiser les résultats.

Le fichier vrt est alors le suivant:

<OGRVRTDataSource>
    <OGRVRTLayer name="BUR Burnot">
        <SrcDataSource>/chemin/varform.shp</SrcDataSource>
        <SrcSQL>SELECT * from varform WHERE "VAR_FORM" = "BUR"</SrcSQL>
    </OGRVRTLayer>
</OGRVRTDataSource>

Et le résultat est:

En pratique, tous les exemples illustrés dans SQL dans OGR peuvent être utilisés.

découper une zone d'un shapefile en fonction d'un polygone;

<OGRVRTDataSource>
    <OGRVRTLayer name="clip">
        <SrcDataSource>/chemin/varform.shp</SrcDataSource>
        <SrcRegion clip="true">POLYGON((170504.41570972 116282.08448781,170576.34282934 113909.22348016,175203.46747074 114069.66549297,175169.19847681 116087.90004614,170504.41570972 116282.08448781))</SrcRegion>
    </OGRVRTLayer>
</OGRVRTDataSource>

Résultat:

il y a aussi moyen d'effectuer des opérations géométriques simples comme l'union entre 2 couches:

Le fichier est alors le suivant (union des BUR entre les 2 couches):

<OGRVRTDataSource>
	<OGRVRTUnionLayer name="union_couche_BUR">
		<OGRVRTLayer name="couche1">
			<SrcDataSource>/chemin/biesm.shp</SrcDataSource>
			<SrcSQL>SELECT * from biesm WHERE "VAR_FORM" = "BUR"</SrcSQL>
		</OGRVRTLayer>
		<OGRVRTLayer name="couche2">
			<SrcDataSource>/chemin/goz.shp</SrcDataSource>
			<SrcSQL>SELECT * from goz WHERE "VAR_FORM" = "BUR"</SrcSQL>
		</OGRVRTLayer>
	</OGRVRTUnionLayer>
</OGRVRTDataSource>

Résultat:

quelques autres fonctions sont aussi disponibles comme la sélection en fonction des surfaces/aires:

<OGRVRTDataSource>
    <OGRVRTLayer name="surfaces">
        <SrcDataSource>/chemin/varform.shp</SrcDataSource>
        <SrcSQL>SELECT * FROM varform  WHERE OGR_GEOM_AREA > 10000000</SrcSQL>
    </OGRVRTLayer>
</OGRVRTDataSource>

Résultat:

SQLite SQL dialect

Ce nouveau format, apparu dans la version 1.10 de GDAL, permet d'utiliser le SQL utilisé par SQLite (SQL As Understood By SQLite), mais aussi toutes les fonctions spatiales de SpatiaLite (SpatiaLite 4.0.0: SQL functions reference list ).

A titre de démonstration, nous allons appliquer plusieurs traitements à une même couche (simplification, zone tampon, enveloppe convexe (Convex hull):

<OGRVRTDataSource>
    <OGRVRTLayer name="BUR">
        <SrcDataSource>/chemin/bur.shp</SrcDataSource>
    </OGRVRTLayer>
    <OGRVRTLayer name="BUR_simplifié">
        <SrcDataSource>/chemin/bur.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT Simplify(geometry,10) from bur</SrcSQL>
    </OGRVRTLayer>
    <OGRVRTLayer name="BUR_tanpon">
        <SrcDataSource>/chemin/bur.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT Buffer(geometry,200) from bur</SrcSQL>
    </OGRVRTLayer>
    <OGRVRTLayer name="BUR_convexhull">
        <SrcDataSource>/chemin/bur.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT ConvexHull(geometry) from bur</SrcSQL>
    </OGRVRTLayer>
</OGRVRTDataSource>

L'ouverture de cette couche vrt dans QGIS conduit au dialogue suivant:

Et en sélectionnant toutes les couches, le résultat est:

Attention, tout n'est pas possible. Ainsi si je veux extraire les limites de polygones (ExteriorRing()), cela ne marchera pas s'il y a plusieurs polygones dans une couche:

<OGRVRTDataSource>
    <OGRVRTLayer name="multipolygones">
        <SrcDataSource>/chemin/bur.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT ExteriorRing(Geometry) from bur</SrcSQL>
    </OGRVRTLayer>
    <OGRVRTLayer name="polygone">
        <SrcDataSource>/chemin/bur1.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT ExteriorRing(Geometry) from bur1</SrcSQL>
    </OGRVRTLayer>
</OGRVRTDataSource>

Seul le deuxième résultat donnera une couche valide. L'autre ne fournira qu'une table non spatiale:

Ceci veut dire que l'on peut parfaitement se servir de ces formats virtuels pour créer des tables non spatiales à sa guise (pour les jointures, etc.).

À quoi cela peut-il bien servir ?

Tout ça, c'est intéressant, me direz-vous, mais à quoi cela peut-il servir en pratique puisque tout peut être fait par un logiciel SIG ?

Cela permet de « bidouiller » ce que l'on veut à partir de couches vectorielles avant traitement dans un logiciel SIG ou pour les traitements sans logiciel SIG:

Avec GDAL/OGR seul:

transformation directe en un autre format:

ogr2ogr -f KML -t_srs EPSG:4326 test.kml test.vrt

création de grilles à partir de points avec divers algorithmes d'interpolation avec la commande gdal_grid (voir par exemple How to generate Concave Hull from regularly spaced point data using GDAL ou LIDAR data to raster file with open source GDAL tool )

gdal_grid -zfield "z" -a invdist:power=2.0:smoothing=1.0 -txe 85000 89000 -tye 894000 890000 *-outsize 400 400 -of GTiff -ot Float64 -l dem test.vrt test.tiff

création de PDF spatiaux (Geospatial PDF) à partir d'OpenStreetMap: voir Making geospatial PDF maps from OSM with GDAL (fichier pdf). Un des fichiers résultants peut être téléchargé à Geospatial_OpenStreetMap_vector_and_raster_map.pdf.

Avec Python:

possibilité de créer ou de manipuler ces fichiers (voir ogr2vrt.py);
comme c'est un format XML, il y a aussi moyen d'utiliser le module ElementTree pour créer un fichier vrt à partir des fichiers shapefiles présents dans un dossier, par exemple;
d'autres traitement sont aussi possibles comme dans Creating a grid from scattered data using inverse of the distance with python (gdal_grid approach).

Avec d'autres programmes:

The power of GDAL virtual formats montre comment utiliser ces fichiers pour servir les données avec MapServer;
Il est aussi possible de les utiliser avec Mapnik, avec FME (FME Supports OGR Virtual Dataset (VRT)), avec les Google Fusion Table (voir Google Fusion Table, OGR et QGIS), etc..