PyQGIS: des géométries, des vecteurs, de l'algèbre vectorielle ou des cosinus directeurs...

dim 02-06-2013 Martin Laloux

SIG OpenSource


	Niveau	Intermédiaire
	Logiciels utilisés	Quantum GIS (QGIS)
	Plateforme	Windows \| Mac \| Linux \| FreeBSD

Traiter des points, des lignes ou des courbes, ça ne vous rappelle rien ? C'est un domaine que vous avez pourtant abordé lors de vos études secondaires dans les cours de géométrie vectorielle et de trigonométrie (du moins je l'espère...).

Quel rapport avec les logiciels SIG ? Ils traitent le même genre de données (on parle bien de couches vectorielles) et donc les mêmes solutions peuvent leur être apportées.

Si je regarde les questions qui sont posées sur les forums géomatiques ou sur gis stackexchange, il est possible de constater que, dans ce domaine, c'est toujours les mêmes questions qui reviennent (car la plupart ne peuvent pas nécessairement être résolues avec les menus ou les boutons du logiciel utilisé ou par un plugin, donc problème...):

comment créer une ligne à partir d'un point et un azimut;
comment obtenir l'azimut d'une ligne;
comment calculer l'angle entre 2 lignes;
comment calculer la distance minimum entre un point et une ligne;
comment obtenir la droite perpendiculaire entre un point et une ligne;
comment créer des points médians ou des points équidistants sur une ligne;
comment découper une ligne en segments;
comment appliquer un changement d'échelle, une translation ou une rotation à un objet;
etc.

Les algorithmes nécessaires pour répondre à toutes ces questions ont pourtant été développés par, notamment, Paul Bourke dès 1998 (paulbourke.net/geometry/) pour les objets 2D et 3D. Il en a même proposé des implémentations en divers langages.

De surcroît, les vecteurs sont utilisés dans bien d'autres domaines et il y a une multitude de tutoriels et de cours (universitaires ou non) qui expliquent, analysent et/ou décortiquent le sujet et apportent des solutions. Il suffit d'aller se servir, mais personne ne semble faire le lien, chacun restant cloisonné dans son domaine...

Comme géologue, j'utilise depuis longtemps des vecteurs dans d'autres spécialités que la géomatique, surtout ceux qui concernent les objets 3D. Mais, malheureusement, QGIS, avec PyQGIS, ne le permet pas encore (voir QGIS, représentation 3D des couches vectorielles (shapefiles dits 3D ou shapefiles avec attributs z) avec les modules Python Matplotlib ou Visvis à partir de la console Python ). C'est pourquoi j'emploie, de préférence, le module Python Shapely.

Après avoir présenté ici les principes, je les appliquerai dans un article suivant pour illustrer quelques applications sur un plan 2D avec PyQGIS à partir des diverses réponses fournies par moi et d'autres sur le ForumSIG et sur gis stackexchange.

Je ne parlerai pas ici des attributs pour me focaliser sur les simples objets géométriques et ne pas créer de difficultés supplémentaires. Je précise aussi que je travaille dans un système de projection cartésien où la géométrie euclidienne a un sens (distances, etc.).

La géométrie avec PyQGIS

Comme déjà soulignés (Python: traitement des couches vectorielles dans une perspective géologique, lecture et enregistrement des couches sous forme de dictionnaires avec le module Fiona), les tutoriels sur PyQGIS ne brillent pas par leur clarté (pour un amateur de Python et d'autres modules géospatiaux), poussant, de plus, les gens a créer tout de suite des plugins, sans bien maîtriser PyQGIS.

Pourtant, le traitement des objets géométriques n'a, en soi, rien de bien compliqué (www.qgis.org/pyqgis-cookbook/geometry.html).

L'objet fondamental est le point, caractérisé par une paire de coordonnées x, y. Une ligne ou un polygone sont composés de plusieurs points:

un point en géométrie PyQGIS se construit avec la commande: QgsPoint(x,y);
la construction d'une ligne est basée sur ces points: QgsGeometry.fromPolyline([QgsPoint(x1,y1),QgsPoint(x2,y2)]));
de même qu'un polygone qui doit être fermé: QgsGeometry.fromPolygon([[QgsPoint(x1,y1),QgsPoint(x2,y2), QgsPoint(x3,y3),QgsPoint(x1,y1)]]);
une polyligne est composée de plusieurs lignes et donc de plusieurs points;
etc.

Les vecteurs

En géométrie euclidienne, les vecteurs ressemblent à ces éléments de géométrie: ils peuvent aussi être définis par une paire de coordonnées x,y,z s'ils décrivent une position dans l'espace. La différence réside dans l'interprétation géométrique et les opérations qui s'y rattachent. Pour simplifier, un vecteur est considéré comme une flèche ayant pour extrémités un point de départ et un point d'arrivée.

Dans le cas des points QGIS (plan en 2 dimensions), ils peuvent être considérés:

soit comme des vecteurs nuls, vecteurs particuliers dont la longueur est nulle, c'est-à-dire dont l'origine est confondue avec l'extrémité;
soit comme des vecteurs normaux avec une origine (point(0,0)) et une extrémité (coordonnées x,y du point).

C'est cette dernière conception qui sera utilisée, car beaucoup plus intéressante comme nous allons le voir. Dans la suite, ne vous étonnez pas du type de coordonnées utilisé (0.5, etc.) mais, toutes les figures ont été réalisées sur QGIS avec le script présenté dans la deuxième partie, crea_mem_layer.py et s'il est appliqué à des coordonnées plus réalistes, le résultat sera moins clair pour les explications (étiquettes, etc.).

points QGIS (QgsPoint(0.5,-1), ...) et vecteurs correspondants:

Vecteur correspondant à la ligne QgsGeometry.fromPolyline([QgsPoint(0.5,-1.0),QgsPoint(1.5,0.5)];

Les coordonnées d'un vecteur sont obtenues par la paire x2-x1, y2-y1, notée entre crochets. Dans le cas du vecteur équivalent à un point QGIS, ce sont donc les mêmes que celles du point (x-0, y-0).

Suivant cette conception, il n'y a donc aucune différence entre un vecteur « point » de QGIS et un vecteur « ligne », tous deux ayant une longueur et une direction.

Définitions

Mathématiquement, un vecteur est un élément d'un espace vectoriel qui répond à un certain nombre de propriétés (fr.wikipedia.org/wiki/Espace_vectoriel#D.C3.A9finitions) et les nombreux cours ou tutoriels qu'on trouve sur Internet. Il est formellement défini par trois composantes:

sa direction
son sens ou son orientation
sa longueur nommée magnitude ou norme

De ce fait, un vecteur peut être déplacé n'importe où dans l'espace, à condition que ses trois composantes, direction, sens et longueur ne changent pas. Ainsi tous les vecteurs de la figure suivante sont égaux, même si les positions des points sont différentes, à la différence des vecteurs géomatiques qui ne sont que des cas particuliers nommés vecteurs de position ou vecteurs de coordonnées (il y a bien d'autres domaines ou la notion de vecteur n'a aucune composante de position comme les vecteurs en physique et autres disciplines, vecteurs de forces, de contraintes, de déformation, etc., voir une très belle illustration de ces concepts en géologie dans les articles de Pierre Thomas sur Planet-Terre, Une bélemnite cisaillée : déformation et contraintes ou En avril, les poissons se font failler):

Vecteur unitaire ou normalisé, vecteurs de position, vecteurs de direction

Un vecteur unitaire ou vecteur normalisé est simplement un vecteur dont la longueur ou la magnitude est égale à 1 (en rouge, vecteur unitaire de la ligne bleue).

Quel est son intérêt ?

il va permettre de nombreux traitements comme la création de points équidistants sur une ligne (simple multiplication du vecteur unitaire).
de plus, comme le même vecteur peut être défini de diverses manières:
- par des coordonnées cartésiennes x,y, (z) (vecteur de position);
- par des coordonnées polaires (d,θ) ou des coordonnées sphériques (d,φ, θ) (vecteur de direction);
- par les cosinus directeurs cosa,cosb, cosc (d.cos(α), d.cos(β),(d.cos(γ)) (vecteur de direction);

Il est facile de constater que dans ce cas:
- il est possible à tout moment de passer des coordonnées cartésiennes à des coordonnées angulaires de direction (attention, sans indication sur le sens ou l'orientation du vecteur);
- le vecteur devient, en coordonnées polaires: x=cos(θ) et y= sin(θ);
- et avec les cosinus directeurs : x = cos(α) et y = cos(β).

Cette particularité est surtout utilisée dans le cas des vecteurs 3D car la direction univoque obtenue peut directement être utilisée, en géologie par exemple, pour tous les traitements d'angles:

En résumé

Ce qu'il faut retenir, c'est qu'à l'inverse des objets purement géométriques, il est possible d'ajouter, de soustraire, de multiplier ou de diviser ces vecteurs et toutes ces opérations ont un sens géométrique simple. Sans le savoir, la plupart d'entre nous effectuent des traitements vectoriels dans nos logiciels SIG mais, de manière non formalisée (voir The mathematics of GIS de Wolfgang Kainz, fichier PDF).

Comme géologue, la meilleure approche du sujet que j'ai eu vient de deux articles proposés dès 1999 par H.L. Vacher dans la revue « Journal of Education Geology » devenue « Journal of Education Geoscience » dans la série Computational geology avec notamment (fichiers PDF):

CG-4, Mapping with vectors, 1999, Journal of Geoscience Education, v. 47, n. 1, p. 64-70;
CG-15, More Mapping with Vectors, 2001, Journal of Geoscience Education, v. 48, n. 3, p. 193-204.

L'application géologique des vecteurs sous forme de cosinus directeurs est amplement illustrée par E. de Kemp pour la modélisation géologique 3D (voir Grass Gis et Paraview: la modélisation géologique 3D pour les démunis, ou les enthousiastes...), R. Groshong dans 3-D Structural Geology (Springer, 2008,) ou R. Allmendinger et N. Cardozo dans Structural Geology Algorithms Vectors and Tensors (Cambridge University Press, 2011) avec, en particulier, l'analyse des tenseurs.

Algèbre vectorielle avec PyQGIS

Qui dit algèbre vectorielle en Python pense à Numpy mais l'utiliser implique de transformer les données en tableaux (array) numpy, ce qui n'est pas le but ici. Il faut plutôt se tourner vers des modules équivalents comme Euclid, pyplane ou les très nombreuses classes Vector ou Vecteur que l'on peut trouver sur Internet comme celles de nliautaud.fr/wiki/articles/vecteurs, code.activestate.com/recipes/52272-vector-a-list-based-vector-class-supporting-elemen/, gist.github.com/mcleonard/5351452, etc., qui n'ont, à priori, aucun lien avec les traitements géomatiques (et oui, tout existe déjà, il suffit de se servir...). Nous n'utiliserons que le module math, présent dans la distribution standard de Python, et le sujet sera abordé sans formalisme mathématique avancé pour ne pas « noyer » le lecteur (comme dans mon cas...).

Formellement pourtant, il faudrait traiter des vecteurs et non des géométries QGIS. Mais, comme une ligne QGIS allant de (x1, y1) à (x2, y2) est représentée par le vecteur [x2 - x1, y2- y1], il est possible de « tricher » et d'assimiler celui-ci à un point virtuel QgsPoint((x2 - x1, y2- y1) pour les calculs, même s'il ne sera jamais figuré.

Longueur d'un vecteur:

La longueur (magnitude, norme) est calculée par une fonction euclidienne classique et, du fait de la définition retenue, un vecteur « point » de QGIS a une longueur (de l'origine (0, 0) au point (x, y)):

def mag(point):
     #magnitude d'un vecteur
     return math.sqrt(point.x()**2 + point.y()**2)

C'est la même chose avec un vecteur « ligne » (la fonction précédente n'est qu'une simplification de dist(point, origine)):

def longueur_theorique(point1,point2):
    #théorique: distance euclidienne entre 2 points
    return math.sqrt((point2.x()-point1.x())**2 + (point2.y()-point1.y())**2)

mais, comme PyQGIS possède une fonction pour le faire (sqrDist, carré de la distance):

def dist(point1,point2):
    #PyQGIS: distance euclidienne entre 2 points#
    return math.sqrt(point1.sqrDist(point2))

Azimut d'un vecteur

L'azimut d'un vecteur se calcule aussi suivant la fonction euclidienne avec ajustements en fonction de l'intervalle voulu (de 0 à 360°):

def azimut_theorique(point1, point2):
    #théorique: calcul de l'azimut entre 2 points,intervalle de 0 à 180°,ici
     angle = math.atan2(point2.x() - point1.x(), point2.y() - point1.y())     
     return math.degrees(angle) if angle &gt; 0 else math.degrees(angle) + 180

Encore une fois PyQGIS possède une fonction pour le faire (azimuth, résultat directement en degrés):

def azimut(point1,point2):
   # PyQGIS: azimut entre entre 2 points (en degrés)
   return point1.azimuth(point2)

d'où l'azimut d'un vecteur « point » ( = azimut(origine, point))

def azimut_pt(point):
    origine = QgsPoint(0,0)
    return origine.azimuth(point)

Addition, soustraction, multiplication

Ce sont des fonctions triviales très faciles à implémenter:

def diff(point2, point1):
     # soustraction de 2 vecteurs
     return QgsPoint(point2.x()-point1.x(), point2.y() - point1.y())

def som(point2,point1):
     # addition entre 2 vecteurs
     return QgsPoint(point2.x()+point1.x(), point2.y()+point1.y())

Il est aussi possible d'effectuer des sommes ou des soustractions avec un scalaire:

def som_k(point, dx,dy):
     return QgsPoint(point.x()+dx, point.y()+dy)

De la même manière, un vecteur peut aussi être multiplié ou divisé par un scalaire (changement d'échelle):

def vecxscal(point,k):
     # multiplication d'un vecteur par un scalaire
     return QgsPoint(point.x()*k, point.y()*k)

PyQGIS a une fonction équivalente (multiply) mais elle traite l'objet géométrique lui-même au lieu de fournir un nouvel objet comme résultat (= vecteur original x k) :

def vecxscal_qgis(point,k):
     # PyQGIS: multiplication d'un vecteur par un scalaire'''
     return QgsPoint.multiply(k)

Norme et normalisation d'un vecteur

Pour rappel, normaliser un vecteur consiste à ramener sa longueur à 1 (vecteur unitaire):

def norm(point1, point2):
    # normalisation d'un vecteur
    point = diff(point2,point1)
    m = mag(point)
    return QgsPoint(point.x()/m, point.y()/m)

et pour un simple vecteur « point » (= norm(origine, point)):

def normpt(point):
    # normalisation d'un vecteur point'''    
    m = mag(point)
    return QgsPoint(point.x()/m, point.y()/m)

Produit scalaire

Le produit scalaire (dot product) entre 2 vecteurs donne, comme son nom l'indique, un scalaire comme résultat (voir Produit scalaire). En pratique, celui-ci est le cosinus de l'angle entre les 2 vecteurs multiplié par les longueurs de ces vecteurs. Parmi ses nombreuses propriétés, il permet donc de calculer l'angle entre les 2 vecteurs et leur orthogonalité: il est égal à 0 si les vecteurs sont perpendiculaires.

def dot_product(point1, point2):
    # produit scalaire de 2 vecteurs 
    return (point1.x()*point2.x() + point1.y()*point2.y())

Calcul de l'angle entre deux vecteurs:

a = QgsPoint(2,3)
b = QgsPoint(5,4)
math.degrees(math.acos(dot_product(a,b)/(mag(a)*mag(b))))
17.650124219930088

Mais avec PyQGIS, ceci peut être fait plus simplement avec:

azimut_pt(b) - azimut_pt(a)
17.65012421993012

Produit vectoriel et déterminant

Le produit vectoriel (cross product) est uniquement défini dans l'espace euclidien tridimensionnel. Le résultat est en effet un nouveau vecteur qui est perpendiculaire au plan contenant les 2 vecteurs (c'est le pôle d'un plan en géologie). Il peut être calculé avec Shapely mais, pas avec PyQGIS (pas d'accès à la coordonnée z).

Avec Shapely:

def cross_product(point1, point2):
      # produit vectoriel avec shapely
      return Point(point1.y*point2.z - point1.z*point2.y, point1.z*point2.x - point1.x*point2.z, point1.x*point2.y - point1.y*point2.x)

Mais, que se passe-t-il si nous considérons 2 vecteurs situés sur un plan (2D, z=0) dans l'espace 3D;

A = Point(2,3,0)
B = Point(5,4,0)
cp = cross_product(A,B)
print (cp.x, cp.y, cp.z)
(0.0, 0.0, -7.0)

Le résultat sera sera bien un vecteur perpendiculaire à ce plan (d'où x= 0 et y = 0) où seule la valeur z a un sens mais, apparemment peu utile avec PyQGIS.

Le même résultat peut être obtenu avec le déterminant de ces vecteurs dans l'espace 2D. Son résultat est un scalaire égal à cette valeur z (je vous laisse le soin de trouver la démonstration mathématique de cette propriété) et cette fonction peut être implémentée dans PyQGIS:

def det(point1, point2):
      # déterminant de 2 vecteurs
      return (point1.x*point2.y) - (point1.y*point2.x)
det(A,B)
-7.0

Et ce déterminant a une signification, voir Déterminant de deux vecteurs dans le plan euclidien dans Déterminant (mathématiques):

il permet de vérifier si 2 lignes sont parallèles (vecteurs colinéaires): résultat = 0;
il permet aussi de calculer l'angle entre les vecteurs;
il permet de calculer l'aire de la surface comprise entre ces vecteurs.

Coordonnées polaires, cosinus directeurs

Transformation vecteur(x, y) vers vecteur (d, θ) en fonction de l'azimut

Il s'agit de la transformation des coordonnées cartésiennes en coordonnées polaires. la manière classique se fait avec la pente: atan(dy/dx):

def car_pol_théorique(point):
    r = mag(point)
    angle = math.degrees(math.atan2(point.y(),point.x()))
    return r, angle

mais comme:

car_pol_theorique(QgsPoint(5,4))
(6.4031242374328485, 38.659808254090088)
90 - azimut_pt(QgsPoint(5,4))
38.659808254090088

il est possible de le calculer en fonction de l'azimut (angle = atan(dx/dy), directement accessible avec PyQGIS:

def car_pol(point):
    return mag(point),90 - azimut_pt(point)

et l'inverse:

def pol_car(dist, angle):
      return QgsPoint(dist * math.cos(math.radians(angle)),dist * math.sin(math.radians(angle)))

vecteur(x, y) vers vecteur(cosa, cosb), cosinus directeurs

Pour rappel, les cosinus directeurs d'un vecteur sont les cosinus des angles entre le vecteur et chaque axe (cosa = d.cos(α), cosb = d.cos(β)).

def cosdir(point):
   # cosinus directeurs d'un vecteur point(x,y)
   cosa = point.x() / mag(point)
   cosb = point.y()/ mag(point)
   return cosa,cosb

Le résultat est un pur vecteur de direction:

cosa,cosb = cosdir(QgsPoint(5,4))
(0.55470019622522915, 0.83205029433784372)
# retour à l'azimut et aux coordonnées polaires définies précédemment
math.degrees(math.atan2(cosa,cosb))
51.340191745909912 # azimut du vecteur point = azimut_pt(QgsPoint(5,4))
math.degrees(math.atan2(cosb/cosa))
38.659808254090088 # angle θ des coordonnées polaires

Il y a aussi moyen de calculer des cosinus directeurs à partir de l'azimut et du pendage en 3D d'une mesure structurale géologique, méthode utilisée, entre autres, par E. de Kemp pour la modélisation géologique 3D mais, comme PyQGIS n'autorise pas la manipulation de la valeur z, il n'est possible de le faire qu'avec l'azimut:

def cosdir_azim(azim):
    az = math.radians(azim)
    cosa = math.sin(az)
    cosb = math.cos(az)
    return cosa,cosb

Le résultat donnera toujours un vecteur unitaire, ce qui facilite grandement la vie pour une multitude de traitements. Le résultat avec des cosinus élimine les petits problèmes rencontrés avec la fonction tangente.

Rotations, translations, cisaillements

Il est possible d'effectuer tous ces traitements avec PyQGIS, il suffit d'aller se servir dans les nombreuses classes Vector ou Vecteur que l'on trouve sur Internet (voir aussi Rotation vectorielle). Cela peut se faire avec les coordonnées cartésiennes ou avec les cosinus directeurs (matrices 2D de translation, de rotation et de cisaillement, sans utiliser Numpy). Voici un exemple de rotation orthogonale et de rotation isométrique après translation de la ligne bleue, avec les nouvelles possibilités d'étiquetage de la version QGIS master:

Mais, depuis que la dernière version du module Shapely dispose de toutes ces fonctions (github.com/Toblerity/Shapely/blob/master/shapely/affinity.py), je l'utilise. Ce sujet étant rarement abordé dans les discussions sur QGIS, je ne le traiterai pas ici (il existe un plugin pour le faire Affine Transformations mais, il est difficile à utiliser si on n'en connait pas bien les principes). Notez que toutes ces opérations permettent de « géoréférencer » une couche vectorielle.

Conclusions

Nous voici armés pour la deuxième partie, l'application tous ces éléments dans PyQGIS: des géométries, des vecteurs, de l'algèbre vectorielle ou des cosinus directeurs, exemples d'applications. Toutes les fonctions ont été regroupées dans le script algèbre_vect_PyQGIS.py.

Tous les traitements ont été effectués sur Mac OS X, avec le script crea_mem_layer.py (origine des figures, hormis celles en 3D) et plusieurs versions de QGIS master.