[QGis 2.x – Spatialite 4.x] Reprojeter des points sur un axe.

Objectif : Suite à la création d’un profil d’après un MNT (voir https://archeomatic.wordpress.com/2013/12/12/qgis-2-0-1-profile-tool-points2one-creer-une-coupe-dapres-un-mnt/) Nous disposons d’un axe de coupe et de mobilier archéologique côté sous forme de points en X,Y et Z. Nous allons « reprojeter » ces points selon l’axe de coupe longitudinal de la fosse. 

 

Reprojection du mobilier archéologique (points) sur un axe de coupe.

Reprojection du mobilier archéologique (points) sur un axe de coupe.

La solution proposée ici consiste a reporter les points sur l’axe de coupe, a calculer leur distance au point 0 (ici le point A) puis a les afficher dans un repère en X, Z (X étant la distance au point 0 et Z, l’altitude du point).

 

Nous allons pour cela utiliser les opérateurs spatiaux SQL de Spatialite (voir la présentation du logiciel dans ce précédent article du blog : http://wp.me/p29w2o-4I ) et notamment st_line_locate_point et st_point.

 

1)      Postulat :

 

Nous disposons dans notre projet Qgis :

–          d’un shapefile de points qui représente le mobilier archéologique avec une valeur Z en attribut.

–          D’un shapefile de ligne qui correspond à une coupe longitudinale. (Note : c’est ce même axe qui a servi a créer un profil d’après le MNT sous-jacent dans le tuto précédent).

 

2)      Créer une base de données Spatialite:

 

Ayant toujours des bugs avec Qspatialite je passe toujours par Spatialite lui-même pour créer une base de données spatialite… rien de plus simple il suffit d’ouvrir l’executable « spatialite_gui.exe) téléchargé ici : https://www.gaia-gis.it/fossil/spatialite_gui/index

 

Puis d’appuyer sur l’icône 2-Connect_create a new DB

Creating a new empty SQlite DB.

 

Le logiciel nous propose alors de créer une nouvelle base de données que nous nommerons REPROJECTION.sqlite

 

 

 

Et voilà c’est fini pour Spatialite, nous pouvons quitter le logiciel et revenir à QGis pour…

 

 

 

3)      Se connecter à la base de données Spatialite et charger nos 2 couches :

 

Dans QGis donc, cliquer sur l’icône 3-Icone-Ajouter-une-couche- , une fenêtre de dialogue s’ouvre alors :

4-Ajouter une couche spatialite

Il suffit d’appuyer sur le bouton [Nouveau] de parcourir pour trouver notre base spatialite nouvellement créée puis d’appuyer sur le bouton [Connecter].

 

Tout le reste des manipulations se feront dorénavant à partir du gestionnaire de Base De Données integré à Qgis 4b-icone Gestionnaire de DB

(dans la barre d’outils « Base de données, la bien nommée !).

 

La fenêtre de dialogue du gestionnaire de Base de données s’ouvre :

5-Gestionnaire-DB---Importe Pour se connecter à notre base de données, il suffit de dérouler le menu spatialite  1 (appuyer sur le petit + ) puis de cliquer sur l’intitulé REPROJECTION.sqlite. Notre base de données est déjà prête et contient une « foultitude » de table par défaut… mais pas encore de données spatiales.

Pour cela nous allons cliquer sur l’icône « Importer une couche ou un fichier » 2

6-Importer une couche vecteur

Il faudra ici :

–          Saisie : choisir notre couche (parmi celle présente dans Qgis) ici Mobilier.shp

–          Table : lui donner un nom : ici SP_Mobilier (le préfixe SP pour Spatialite qui permet de la differencier de l’originale)

–          Clé primaire : cocher et écrire gid

–          Colonne de géométrie : cocher et écrire geom

–          Créer des géométries simples au lieu de multipartities : cocher

–          Créer un index spatial : cocher (ça ne peut pas faire de mal !)

Appuyer sur le bouton [OK]

 

Et voilà une belle couche importée. Pour le vérifier, cliquer sur l’intitulé de la base de données puis rafraichir (icône entouré de rouge ci-dessous). On voit bien que notre couche de point SP_Mobilier a été importée :

7-couche importée

On fera de même avec notre coupe Coupe.shp que l’on appellera SP_Coupe.

8_Apercu-table

Une fois les couches importées on peut très simplement visualiser nos couche dans l’arborescence 1 et la table attributaire de nos couches en cliquant sur l’onglet « Table » 2.

On voit par exemple ci-dessus que nos points ont comme champs attributaires : gid (la clé primaire), geom (la colonne de geométrie), x,y et z (no comment), type et code (description du mobilier archéo) et fait (le numéro de fait archéologique qui contient le mobilier).

 

4)      La reprojection proprement dite…

 

… va se passer en 2 étapes …

Note :

– Boouuh avec PostGis on peut le faire en une fois !!

–   oui oui mathias je sais –et merci au passage pour ton aide indispensable- mais là pas besoin d’instalation client-serveur compliquée.. et toc !

et le tout dans la fenêtre SQL à laquelle on accède via l’icône 9-Icone fenetre SQL .

 

4a) Créer une table de points avec un « xrelatif » :

10-Creer-table-de-points

Nous donnons donc ici un nom 1 à notre première requête « Créer_Table_Point » puis après l’avoir tapée 3 le bouton [Stocker] 2 permettra de la réutiliser.

Avant d’appuyer sur le bouton [Executer] 4 et d’avoir vérifier que la requête SQL avait été traitée convenablement 5 attardons nous sur la requête :

 

create table T301 as

on crée une table T301 (T pour table et 301 car c’est le nom de la fosse) avec

Select p.gid, p. »z », p. »type », p. »code », p. »fait »,

les champs gid, z, type, code et fait…(notez les noms de champs avec les guillemets)

st_line_locate_point(c.geom, p.geom)*st_length(c.geom) as xrelatif

… et un champ appellé xrelatif qui correspond à la distance au point 0 de la coup, de chaque point préalablement reporté sur l’axe de coupe

from « SP_Coupe » c, « SP_Mobilier » p

et ce, à partir des couches SP_Coupe (c de son petit nom) et SP_Mobilier (p pour point de son petit nom)

where p. »fait » like ‘301’ and  c. »gid » = 1

avec comme condition que les points appartiennent au fait 301 et que la clé primaire de la coupe (gid) soit égale à 1.(notez la valeur texte encadrée d’apostrophes)

order by xrelatif ;

et trions tout par ordre croissant à partir du champs xrelatif !

 

Vous avez tout compris ? Vous êtes donc autorisé a [Executer] votre requête via le bouton idoine …choisir le SCR… (en passant si la requête a fonctionné appuyer sur le bouton [Stocker] pour la garder en mémoire !)…

…et a regarder la jolie table que vous venez de créer… après avoir actualisé votre base de donnée et cliqué sur l’onglet « Table » bien sûr !

11-Table 301

 

Nous avons donc bien une table avec les champs voulus et un nouveau champ : « xrelatif » qui correspond pour chaque point reporté sur l’axe a sa distance par rapport au départ de l’axe (en mètre bien sûr, mon projet QGis et mes couches sont en RGF93/ CC47)

Ne reste plus qu’à…

 4b) Créer une couche de point en Xrelatif, Z :

12-creer-couche-de-points

Rien de plus simple (mais si vous allez voir…)

 

1 Donner un nom à notre requête, ici « Créer_Points_Projetés »

2 Ne pas appuyer sur le bouton [Stocker]… attendre d’avoir écrit la requête et vérifié qu »elle fonctionne proprement 😉

3 Ecrire la requête :

 

Par laquelle on demande

select st_point(xrelatif, z)

de créer les points avec comme abscisse les valeurs du champ xrelatif et en ordonnée les valeurs du champ z (attention pas de create table ici, juste un select)

as geom, gid, « z », « type », « code » , « fait », « xrelatif »

et d’y ajouter comme attributs, les champs geom, gid, z, type, code, fait et xrelatif

from T301

et ce, à partir de la table T301, que l’on vient de créer avec la requête précédente

order by xrelatif

et puis, tant qu’à faire de trier à partir du champ xrelatif

 

On [Exécute]

5 On regarde la table résultante…

6 On coche « Charger en tant que nouvelle couche »

7 choisir la clé primaire : gid

8 identifier la colonne qui contient la geométrie : geom

9 donner un nom à la couche, ici « Reprojection_Points »…

10 et on [Charge !]

REPOJECTION_FINAL

Conclusion :

Alors, c’est-ti-pas-beau ? avec le profil en plus…

Cette méthode embarquée (pas de système client serveur) permet simplement de faire ce dont je rêvais depuis pas mal de temps… reprojeter des points sur un axe afin d’avoir une vision « en coupe » des altitudes de différents objets archéologiques. On peut évidemment choisir l’axe de reprojection à l’envie mais surtout cette méthode totalement intégrée au logiciel d’information géographique permet évidemment de garder tous les attributs des points d’origine et donc de les représenter avec la même symbologie en coupe et en plan.

Personnellement j’y vois déjà de multiples usages à des échelles différentes et à des fins différentes.. oui oui géomorphos de tous poils et responsables d’opérations des profondes vallées, c’est à vous que je parle !

En tout cas si ce tuto vous a servi, dites-m’en plus !

 

Epilogue :

Vous avez bien compris la combine qui se rapproche de celle utilisée dans le tuto précédent pour créer un profil selon un axe d’après un MNT, on triche puisque l’on crée des points dans un système Xrelatif, Z au lieu de X,Y.

Un des plus gros problème qui en résulte c’est que le composeur de Qgis refuse d’afficher les objets trop proches de 0,0 (Bug ?), pour palier a ce détail technique il suffit de changer la 3eme ligne de la première requête ainsi : (oui bah vous avez qu’à suivre un peu !)

st_line_locate_point(c.geom, p.geom)*st_length(c.geom)+1000000 as xrelatif

 

Publicités
Tagué , , , , ,

15 réflexions sur “[QGis 2.x – Spatialite 4.x] Reprojeter des points sur un axe.

  1. Zébulon dit :

    Pas mal !

  2. Emeline dit :

    Je viens de tester ton tutoriel pour reprojeter des points sur un axe, c’est super !!
    Merci beaucoup

  3. Emeline dit :

    Oui oui pour de l’archéo. J’ai relevé une coupe au tachéo en visée direct avec laser.
    Et j’ai utilisé ton tuto pour reprojeter mes points en coupe !!

  4. rozenn dit :

    Super le tuto!! J’ai réussi a balancé tout le shape de points (silex) sur ma ligne!! Par contre, peut-on mettre une contrainte en ne sélectionnant que les points qui sont à une certaine distance de la ligne (par ex 50 cm de chaque côté?)..
    Et puis question bête, comment as-tu fait pour la présentation finale (le plan en haut et la coupe en dessous?), c’est dans le composeur d’impression?
    Merci

    • archeomatic dit :

      Salut Rozenn !

      I- Sélectionner les points à 50cm de la coupe:

      1) avec Qgis tu sais le faire 😉
      a) créer un Buffer (tampon) avec une distance fixe de 0.5m (Vecteur>Outils de Geotraitement>Tampon(s)
      b)Requête Spatiale (Vecteur>Requête spatiale) de Points/A l’intérieur/Buffer
      c) enregistrer la sélection comme une couche

      2) Direct dans Spatialite grâce à la fonction PtDistWithin (équivalent de ST_DWithin dans PostGis)tu peux fixer la condition de distance entre 2 géométries:
      dans notre exemple ça donne:
      where p. »fait » like ‘301’ and c. »gid » = 1 and PtDistWithin(c.geom, p.geom, 0.5)
      ou je demande comme condition que la distance entre la coupe (c.geom) et mes points (p.geom) soit égale à 0.5m !!!

      c’est pas la classe ca ?!

      II- Mettre le plan et le profil sur un même document:

      Cela se fait effectivement via le composeur d’impression on ajoute 2 fois une carte avec une subtilité avant quand même : il faut ajouter 1 000 000 aux coordonnées des points projetés pour pouvoir les faire apparaitre dans le composeur sinon il ne veut rien entendre..

      En attendant des nouvelles : Biz
      A+

  5. rozenn dit :

    SUPER!! Trop balèze le tuto!! MERCI encore…

  6. Sof dit :

    Merci Sly 🙂

  7. Lotin dit :

    Bon je suis une véritable andouille une fois sur Qgis, je m’en sors vraiment pas bien et du coup je n’y fait que des choses basiques.
    Je cherche un moyen de projeter sur les axes de coupes des points dont j’ai besoin (un foyer à galets à chauffés pour aiguiller) afin de voir si les altis inf dessinent un creusement qui nous a totalement échappé à la fouille (suis-je assez clair ?). J’arrive à le faire basiquement mais en trichant, je projette sur l’axe des X les Z puis sur l’axe des Y les Z. Seulement, cela ne correspond pas vraiment à mes axes de coupe, plutôt à des axes nord/sud et ouest-est (logique). Je ne sais même pas par quel tuto commencer pour cela. Des pistes ?

    • archeomatic dit :

      Salut lotin,
      D’après ce que tu le dit ce billet de blog est sensé répondre à ta question.. Effectivement il mêle QGIS et spatialite ce qui relève d’un « niveau utilisateur avancé » mais normalement si tu suis le billet à la lettre tu devrais y arriver..
      Au pire si tu n’est pas pressé tu peux m’envoyer tes données (nuage de points et axes de coupes sur lesquels tu veux ré projeter tes galets)..
      Ou vois avec un topo sous Autocad ou un sigiste de ton entourage.

Laisser un commentaire

Entrez vos coordonnées ci-dessous ou cliquez sur une icône pour vous connecter:

Logo WordPress.com

Vous commentez à l'aide de votre compte WordPress.com. Déconnexion / Changer )

Image Twitter

Vous commentez à l'aide de votre compte Twitter. Déconnexion / Changer )

Photo Facebook

Vous commentez à l'aide de votre compte Facebook. Déconnexion / Changer )

Photo Google+

Vous commentez à l'aide de votre compte Google+. Déconnexion / Changer )

Connexion à %s

%d blogueurs aiment cette page :