[SpatiaLite] corruption de l'index spatial R*Tree : soyons prudents !

Aujourd'hui, il est question de corruption d'index R*Tree dans vos BDD SpatiaLite. En effet, bien que très utile pour optimiser les requêtes, il nécessite quelques précautions d'usage.

Ce post est s'inspire du document suivant : http://www.gaia-gis.it/spatialite-2.4.0-5/SpatialIndex-Update.pdf (écrit par le développeur de SpatiaLite)

Voici comment réussir à corrompre un index Spatial R*Tree:

D'abord, quelques rappels pour mieux comprendre le fonctionnement de R*Tree:

• Tout Index Spatial R*Tree est simplement une table (virtual table)
• Le moteur interne de SQLite ne gère pas de façon native les relations liant l'index à la colonne géométrique de la table correspondante.
• Spatialite dispose d'un système de TRIGGERS afin d'assurer la synchronisation de l'index avec la colonne géométrique de la table correspondante.
• Le lien entre l'index et la colonne géométrique repose sur les valeurs de ROWID
• Pour une table SQLite donnée, chaque ligne est identifiée de façon unique par son ROWID
• Si la table contient une clef primaire (PRIMARY KEY), alors le ROWID est totalement lié aux valeurs de la clé primaire.
• Si la table ne contient pas de clef primaire, alors le ROWID correspond simplement au numéro de ligne et n'est lié avec la table que de façon relative

Dans quels cas l'index peut-il être corrompu ?

Considérons une table SANS clé primaire, ayant un index spatial R*Tree. Maintenant, supprimons une ligne de cette table, et compactons la base de donnée (VACUUM): vous venez de corrompre votre index spatial...

• La commande VACUUM compacte la base de donnée et par conséquent réassigne les ROWID des tables
• Pendant un VACUUM, les triggers sont désactivés ( fonctionnement interne de SQLite). Ainsi, la synchronisation avec l'Index Spatial n'est pas possible.
• On se retrouve dès lors avec une table ayant des ROWID modifiés et un index Spatial non modifié: le lien entre la table (colonne géométrique) et l'index spatial est rompu !
• Les requêtes basées sur l'index spatial ne sont dès lors plus valides....

La seule solution a ce problème: Toujours affecter une clé primaire à vos tables avant de mettre en place un index spatial R*Tree !

• La clé primaire va permettre s'assurer que les ROWID ne soient pas réassignés lors du VACCUUM et donc garantir la synchronisation avec l'index spatial R*Tree

BILAN: Méthodologie à suivre pour créer de manière sûre une table avec clé primaire et index spatial:

ex: matable(pk_uid,nom,geometry)

#Création de la table et mise en place de la clé primaire
CREATE TABLE ma table
( pk_uid INTEGER PRIMARY KEY , nom TEXT NOT NULL )

#Création de la colonne géométrique
SELECT addGeometryColumn( "matable" , "geometry", srid , type , dimension )

#Insertion des données
INSERT INTO matable ( pk_uid,nom,geometry) SELECT ..... FROM .....

#Création de l'Index Spatial
SELECT createSpatialIndex("matable", "geometry")

En procédant ainsi, vous êtes sûr de ne jamais corrompre votre index spatial.

[SpatiaLite] Mettre à jour une table avec jointure (UPDATE SET JOIN)

SpatiaLite, encore et toujours !
Pour les bases, je rappelle que c'est ici que ça se passe:
http://www.gaia-gis.it/spatialite-2.4.0-4/spatialite-cookbook-fr/

Aujourd'hui, nous allons apprendre comment mettre à jour une table en utilisant une relation (jointure). IL ne s'agit pas d'une syntaxe spécifique à SpatiaLite, mais à SQLite en général.

Données:
-parcelles (id, surf_lib, geometrie) -type:polygon-
-parcelles_maj(id, surf_lib)

Mise en situation:
Vous gerez un domaine agricole, composé de plusieurs parcelles. Vous disposez d'une table parcelles datant de l'an dernier, et recensant les surfaces libres (disponibles) sur chacune des parcelles. Vous decidez d'actualiser cette table en envoyant un stagiaire réaliser une verification de certaines parcelles. Il reviens vous voir avec une table parcelles_maj . Comment mettre à jour votre table parcelle à partir des informations présentes dans la table du stagiaire ?

Solution:
UPDATE parcelles
SET surf_lib = ( SELECT a.surf_lib FROM parcelles_maj AS a WHERE a.id=parcelles.id)

Quelques explications:
SQLite n'autorise pas l'utilisation de la clause JOIN dans des reqêtes de type UPDATE...SET...
La syntaxe ci dessus en revanche est acceptée et marche correctement.

Une fois le principe compris, il est facile de mettre à jour une table en fonction de plusieurs tables, autant avec des relations attributaires que spatiales.

En cas de soucis, n'hésitez pas à poser des questions,

[SpatiaLite] Soustraction vectorielle

Aujourd'hui nous allons réaliser une soustraction vectorielle avec SpatiaLite

Données:
- table1 (id,geometry) - type: polygone -
- table2 (id,geometry) + Index Spatial R*Tree - type: polygone -

Objectif: Soustraction vectorielle simple: table3=table1-table2
(l'ordre est important: dans notre cas, on va soustraire à la table 2 à la table)

Mise en situation:
Imaginons que:
- Notre table1 représente des territoires à explorer
- Notre table 2 représente des zones inaccessibles (fossés, lacs, etc...)
Notre objectif serait ici de determiner parmis les territoires à explorer (table1), les zones accessibles.

Solution basique:
CREATE TABLE table3 AS SELECT
t1.id AS id1
st_multi( st_difference( t1.geometry,t2.geometry) ) AS geometry
FROM table1 AS t1
JOIN table2 AS t2 ON (
GeometryType( st_multi(st_difference( t1.geometry,t2.geometry) ) ) = "MULTIPOLYGON")

Quelques explications:
Ici, on créée la table3 contenant le resultat de la soustraction vectorielle via la fonction st_difference(). L'utilisation de st_multi() force le type de géométrie créée en multi: ceci est recommandé pour la validation ulterieure de la colonne géométrique via la fonction RecoverGeometryColumn(). De plus, le critère de jointure GeomtryType()="MULTIPOLYGON" permet de ne garder que les résultats de soustraction de type MULTIPOLYGON, écartant ainsi les géométries NULLes (objets de table1 entièrement contenu dans objets de table2) et les lignes (POLYLINESTRING).
Ainsi, la table3 contiendra uniquement des objets de type MULTIPOLYGON. La fonction RecoverGeometryColumn() permet ainsi d'authentifier correctement la colonne geometry.

Cette requête basique marche, mais son fonctionnement n'est pas optimal et les temps de calcul risquent d'être très longs pour des gros jeux de donnée. Essayons de l'optimiser un peu ...

Solution optimisée:

CREATE TABLE table3 AS
SELECT
t1.id AS id1
st_multi( st_difference( t1.geometry,t2.geometry) ) AS geometry
FROM table1 AS t1
JOIN table2 AS t2 ON (
st_intersects( t1.geometry,t2.geometry)
AND
GeometryType( st_multi(st_difference( t1.geometry,t2.geometry) ) ) = "MULTIPOLYGON"
AND
t2.rowid IN (
SELECT pkid FROM idx_table2_geometry
WHERE pkid MATCH RTreeIntersects(
MBRminX(t1.geometry),
MBRmaxY(t1.geometry),
MBRminX(t1.geometry),
MBRmaxY(t1.geometry)
))
)
UNION
SELECT
t1.id AS id1,
st_multi(t1.geometry) AS geometry
FROM table1 AS t1, table2 AS t2
WHERE
NOT st_intersects( t1.geometry,t2.geometry)

Quelques explications:

L'astuce consiste à filtrer au préalable les objets s'intersectant, avant de réaliser le géotraitement (soustraction vectorielle) uniquement sur ces objets:

Dans un premier temps, on selectionne tous les objets de table1 qui intersectent les objets de table2 via la fonction st_intersects() (l'utilisation du R*Tree en rouge). Pour chacun de ces objets, on réalise une soustraction vectorielle via la fonction spatiale st_difference(obj1,obj2).

Dans un second temps, on ajoute au résultat les objets de table1 qui n'intersectent pas ceux de table2. La clause UNION va permettre de fusionner les résultats des deux SELECT en un seul résultat. Ces objets ne sont pas modifiés.

A vous de jouer !

N'hesitez pas à laisser des commentaires en cas de soucis, je tâcherais de vous aider du mieux que je peux.

GPS garmin etrex venture HC /QGIS et Ubuntu 10.04

Afin de préparer une campagne d'échantillonnage terrain, j'ai eu a utiliser un GPS GARMIN etrex venture hc, sous ubuntu 10.04 équipé de QGIS 1.6. Suite aux nombreux problèmes rencontrés, je propose la solution que j'ai trouvé:

Tout d'abord, vous trouverez ici le forum qui m'a le plus aidé:

http://www.gpspassion.com/forumsen/topic.asp?TOPIC_ID=121724

Installation des logiciels:

Dans un terminal tapez:

>sudo apt-get install gpsbabel

>sudo apt-get install gpsd

>sudo apt-get install Gebabbel

>sudo apt-get install viking

Afin de permettre la connection avec le GPS GARMIN:

>sudo gedit /etc/udev/rules.d/51-garmin.rules

Dans la fenêtre qui s'ouvre, tapez:

SYSFS{idVendor}=="091e", SYSFS{idProduct}=="0003", MODE="666"

Enregistrez et redémarrez

Communication avec le GPS

Utiliser viking pour communiquer avec le GPS: http://doc.ubuntu-fr.org/viking

Transformer SHP en GPX

Utiliser QGIS pour transformer une couche SHP (ShapeFile) en GPX (format d'échange GPX, utilisé par viking entre autre):

>charger couche shp dans qgis

>clic-droit sur la couche shp > sauvegarder sous ...

> choisir le format GPX

> sélectionner le système de projection WGS84 (EPSG: 4236)

Uploader votre GPX dans votre GPS GARMIN ETREx VENTURE HC

>ouvrir viking

>créer un "layer" GPS

>sélectionner le port USB dans "data preference" dans le fenetre qui s'ouvre

>menu>ouvrir> ouvrez votre fichier GPX: un nouveau "layer" apparaît

>sélectionnez les éléments à importer (waypoints, routes, tracks) et par glissé-déposé, ajoutez le dans le sous layer "émission " du layer GPS.

>clic droit sur layer GPS-> Emission vers le GPS.

En espérant vous avoir aidé !

[SpatiaLite] tirage aléatoire dans une base de donnée

Pour apprendre les bases de SpatiaLite, direction le cookbook en Francais: ici

Données:
- batiments (id, surface_sol, geometry) - type: polygon - + Index Spatial Rtree
- riviere (id, nomn geometry) - type: ligne -

Objectif: Tirer aléatoirement 5 batiments ayant une surface au sol > 100m² et situés à moins de 2500m du cours d'eau

Solution:
SELECT
a.*
FROM batiments as a
JOIN riviere as b ON (
ST_Distance(a.geometry,b.geometry)<=2500
AND
a.rowid IN (
SELECT pkid from idx_batiments_geometry
WHERE pkid MATCH RTreeintersects(
MBRminX(b.geometry)-2500,
MBRminY(b.geometry)-2500,
MBRmaxX(b.geometry)+2500,
MBRmaxY(b.geometry)+2500,
))
)
WHERE
a.surface>=100
ORDER BY random()
LIMIT 5

Quelques explications:
Requête SELECT ... FROM ... WHERE ... classique permettant de selectionner les bâtiments de plus de 100 m². La clause JOIN...ON... permet de joindre la table riviere, en filtrant uniquement les bâtiments situés à moins de 2500m d'un cours d'eau via st_distance(). Le code SQL en rouge , facultatif, utilise l'index spatial R*Tree de la table batiments afin de préfiltrer les batiments intersectant le MBR des cours d'eau, élargis de 2500 m, et ainsi améliorer la performance de la requête.
Le point le plus important ici réside dans l'utilisation de la fonction ORDER BY random() : ceci va permettre d'afficher les lignes du résultat dans un ordre aléatoire.
La clause LIMIT 5 permet de ne selectionner que les 5 premiers résultats.


L'association de GROUP BY random() et LIMIT x va permettre de simuler un tirage aléatoire de x lignes d'une table. Ceci est particulièrement intéréssant pour préparer des échantillonnages.
L'un des principaux avantages de réaliser le tirage aléatoire directement par SQL, est de pouvoir specifier des critères sur la population initiale.

Encore plus fort ....

Tentons maintenant de réaliser en une seule étape un échantillonnage stratifié.
Parmis les bâtiments situés à moins de 2500m d'un cours d'eau, je veux échantillonner:
- 5 batiments au hasard parmis les bâtiments de plus de 100 m² (>=100m²)
- 2 batiments au hasard parmis les bâtiments de moins de 100 m² (<100m²)

Solution:
SELECT id FROM
( ...requete précédente... )
UNION
SELECT id FROM
( ...requete précédente, en remplaçant >=100 par <100 et LIMIT 5 by LIMIT 2 ... )

La clause UNION va permettre de fusionner les résultats des deux requêtes. Ainsi, votre résultat se composera de l'id de l'ensemble des batiments à échantillonner.

tutoriel spatialite: pourcentage d'intersection

Objectif: determiner le pourcentage de recouvrement entre deux couches SIG

Données:
- table1 (geometry, id)
- table2 (geometry,id) + Index Spatial RTree

Méthode: on va calculer pour chaque objet de la table 1, l'aire (ou le pourcentage) d'intersection avec les objets de la table 2

Solution:
SELECT
table1.id as id1,
table2.id as id2
st_area( st_intersection(table1.geometry,table2.geometry) ) as "aire de recouvrement",
st_area( st_intersection(table1.geometry,table2.geometry) )/ st_area ( table1.geometry)
as "pourcentage de recouvrement"
FROM table1,table2
WHERE
st_intersects(table1.geometry,table2.geometry)
AND
table2.rowid IN (
SELECT pkid FROM idx_table2_geometry
WHERE pkid MATCH RTreeintersects(
MBRminX(table1.geometry),
MBRminY(table1.geometry),
MBRmaxX(table1.geometry),
MBRmaxY(table1.geometry),
)))

Quelques explications:
Clause SELECT ... FROM: on selectionne les colonnes voulues dans les tables.
- st_area(st_intersection()) calcule l'aire de l'intersection
Clause WHERE:
- jointure spatiale: la fonction st_intersects() permet de selectionner les objets des deux tables se chevauchants.
- RTree: Pour chaque objet de la table1, on préfiltre les objets de la table2 intersectant le RTree. Ceci permet d'accelerer grandement la requete

Facile non ?
En revanche, il est plus difficile de selectionner directement par SQL, pour chaque objet de la table 1, l'objet de la table 2 correspondant au recouvrement maximal (utile pour passer les valeurs d'une table dans une autre en fonction du taux de recouvrement). En effet, la clause GROUP BY id1 couplée à MAX(%intersecteion) va ici nous retourner des lignes composées de la façon suivante:
1) Une ligne par objet de la table 1
2) Pourcentage max de recouvrement entre cet objet et les objets de la table 2
3) un objet de la table2 "au hasard" (ne correspondant pas au taux de recouvrement maximal)

Objectif 2: Determiner pour chaque objet de la table 1, l'objet de la table 2 avec lequel le pourcentage de recouvrement est maximal

Solution:

Etape1:
Creer une vue à partir de la requete précédente ( vue = intermédiaire entre table et requete):
CREATE VIEW vue1 AS
SELECT ....requete précédente....

La vue1 contient alors la liste des pourcentages de recouvrement pour chaque couple d'objets (table1 et table2)

Etape2:
Creer une nouvelle vue contenant, pour chaque objet de la table1, la valeur max du pourcentage de recouvrement:
CREATE VIEW vue2 AS
SELECT
vue1.id1 as id1,
MAX( vue1."pourcentage de recouvrement" ) as "pourcentage de recouvrement max"
FROM vue1

Etape3:
Enfin, en utilisant les deux vues précédentes, il devient possible d'extraire de la vue1, les lignes correspondants au % max de recouvrement pour chacun des objets de la table1:
SELECT
*
FROM vue1
JOIN vue2 ON (
vue2.id1=vue1.id1
AND
vue2."pourcentage de recouvrement max"=vue1."pourcentage de recouvrement"
)

Et voila, le travail est terminé ! A priori cela peut paraître compliqué, mais en décortiquant un peu les requetes, on s'apperçoit que c'est à la portée de tout le monde.
De plus, répondre à ce genre de question via des logiciels SIG classiques (MAPINFO, QGIS, ...) s'avère encore plus complexe à mettre en oeuvre ...

Tutoriel SpatiaLite: distance minimale entre les objets de deux couches SIG

Pour ceux qui ne connaissent pas SpatiaLite, c'est ici que ça se passe:
http://www.gaia-gis.it/spatialite-2.4.0-4/spatialite-cookbook-fr/index.html

Données de départ:
- table1 (id,nom,geometry)
- table2 (id,nom,geometry) + Index Spatial R*Tree

Objectif: Pour chaque objet de la table1 , determiner la distance de l'objet le plus proche dans la table2 (dans la limite de 1Km, soit 1000m)

Contrainte: De manière intuitive, la solution serait toute simple:
SELECT
table1.nom as nom1,
table2.nom as nom2,
min(st_distance(table1.geometry,table2.geometry)) as distance
FROM table 1
JOIN table2 ON ( ....dist...)
GROUP BY 1
A priori, cette requete retourne effectivement pour chaque objet de la table1, la distance la plus courte le séparant des objets de la table2. Mais en regardant de plus près, pour chaque ligne l'objet de la table2 retourné par la requète correspond non pas à l'objet le plus proche, mais au premier objet trouvé... ceci est dû au fonctionnement même de GROUP BY. Afin de determiner les couples d'objets les plus proches, il faut ruser en passant par plusieurs étapes.

Solution: en plusieurs étapes

Etape 1: calcul des distances entre tous les couples possibles.
CREATE VIEW distances AS
SELECT
table1.nom as nom1,
table2.nom as nom2,
st_distance( table1.geometry,table2.geometry ) as distance
FROM table1
JOIN table2 ON (
st_distance( table1.geometry,table2.geometry )<=1000
AND table2.rowid IN (
SELECT pkid FROM idx_table2_geometry
WHERE pkid MATCH RTreeintersects(
MBRminX(table1.geometry)-1000,
MBRminY(table1.geometry)-1000,
MBRmaxX(table1.geometry)+1000,
MBRmaxY(table1.geometry)+1000
)))
GROUP BY 1

Analyse de la première étape:
Création d'une VUE ( que l'on pourra réutiliser dans les prochaines étapes - intermédiaire entre une table et une requete )
Requete du type SELECT.....FROM.... classique.
La fonction spatiale st_distance() permet de calculer la distance minimale entre chaque couple.
La sous-requete SELECT permet d'exploiter l'Index Spatial, en préfiltrant pour chaque objet de la table1, les objets de la table2 dont le MBR élargi de 1km intersecte le R*tree.

Etape 2: calcul de la distance minimale séparant chaque objet de la table 1 des objets de la table2

CREATE VIEW distances_min AS
SELECT
nom1 as nom1,
min( distance ) as "distance min"
FROM distances
GROUP BY 1

Analyse de la deuxième étape:
On crée une nouvelle vue à partir de la vue précédente (distances).
min(distance) et GROUP BY 1 vont permettre de selectionner pour chaque objet de la table1, la distance minimale le séparant des objets de la table2.

Etape 3: A chaque objet de la table1, nous associons l'objet de la table2 le plus proche

SELECT
distances.nom1 as nom1,
distances.nom2 as nom2,
distances.distance as distance
FROM distances
JOIN distances_min ON (
distances.nom1=distances_min.nom1
AND
distances.distance=distances_min."distance min"
)

Analyse de la dernière étape étape:
On va selectionner ( SELECT ... FROM ) les lignes de la VUE Distances.
Afin de récuperer, pour chaque objet de la table1, uniquement l'objet de la table2 le plus proche, on va JOINdre la vue distances_min.

Et voila, vous avez maintenant, pour chaque objet de la table 1:
- L'objet de la table 2 le plus proche
- La distance entre ces deux objets