Vous appréciez mon travail ?
Je serais ravi de prendre un café !

Vous prenez du plaisir à lire mes articles ? Vous apprenez de nouvelles choses ? Je serais ravis que vous supportiez mon travail avec une petite participation

1 café Merci, vous financez ma dose quotidienne de théïne (oui, en vrai je ne bois pas de café).
5 cafés Génial, ça couvre mes frais de serveur mensuels.
10 cafés Fantastique, avec ça je peux investir dans du matériel et approfondir mes connaissances.
BazinGa's - Tips & tuto IT

Soyons précis avec PostGIS

Vous êtes vous déjà demandé quelle était la précision de PostGIS ? Moi oui… Alors étudions cela ensemble.

La documentation de PostGIS n’est pas très claire sur la question de la précision mais il faut savoir que toutes les coordonnées sont stockées sous forme d’entiers à virgule flottante en double précision. La précision de PostGIS est donc liée à ce format de stockage des numériques.

Virgule flottante et double précision

Ce paragraphe est une simplification de la norme IEEE 754 pour vous aidez à comprendre le fonctionnement des nombres à virgules flottantes en double precision. Des imprécisions sont donc volontairement présentes.

La virgule flottante est une méthode de représentation d’un nombre en informatique qui s’apparente à la notation scientifique (par exemple 1,234 * 10123). Elle repose ainsi sur les éléments suivants :

  • un signe : + ou -.
  • Une mantisse : une valeur binaire qui définie une fraction entre 0 et 1.
  • Une base : un nombre entier (par exemple 2 ou 10).
  • Un exposant : un nombre entier.

On obtient la valeur finale grâce à l’opération suivante :

x = signe * 1,mantisse * baseexposant

Quand on parle de double précision, on fait alors référence à un format particulier de nombre à virgule flottante :

  • L’encodage se fait sur 64 bits :
    • 1 bit pour le signe (0 = + et 1 = -)
    • 11 bits pour l’exposant
    • 53 bits pour la mantisse dont :
      • 1 bit pour le décalage (voir le paragraphe sur la dénormalisation)
      • 52 bits pour la valeur
  • La base est 2 car le stockage se fait sous forme binaire.
  • L’exposant est décalé, par convention on ajoute 1023 a sa valeur réelle.
  • Le calcul d’un nombre en double précision est alors le suivant :
    x = -1signe * (1 + mantisse) * 2exposant-1023
  • Dans le cas d’un nombre dénormalisé (voir le paragraphe dédié), le calcul est le suivant :
    x = -1signe * (0 + mantisse) * 21-1023 = -1signe * 0,mantisse * 2-1022

Exposant

L’exposant est codé sur 11 bits ce qui permet les valeurs suivantes :

  • 0 : 00000000000 : valeur minimum
    Valeur spécifique réservée pour représenter le 0 ou un nombre dénormalisé (voir plus bas).
  • 1 : 00000000001
  • 2046 : 11111111110
  • 2047 : 11111111111 : valeur maximum
    Valeur spécifique réservée pour représenté les valeurs infini ou NaN (Not A Number : pas un nombre)

Cette valeur d’exposant est non signée. Afin de permettre le stockage d’exposants négatifs, un décalage est introduit. Dans le cas du double précision, ce décalage est de 1023.

Ainsi, avec un exposant décalé de 1023, la plage d’exposants utilisables est la suivante :

  • De 0 - 1023 = -1022
  • À 2046 - 1023 = 1023

Mantisse

La mantisse est codée sur 52 bits et représente une fraction comprise entre 0 et 1 sous la forme d’un nombre binaire. La mantisse peut ainsi avoir 252 (4 503 599 627 370 496 = 4,503599627370496 × 1015) valeurs différentes.

Le calcul de la valeur décimale de la mantisse à partir de la valeur binaire se fait comme suivant :

Valeur binaire : 10100...00101 (la partie tronquée n’est composée que de zéros).

La valeur décimale est égale à la somme des puissances de deux multipliées par la valeur du bit correspondant :
(1*2-1) + (0*2-2) + (1*2-3) + (0*2-4) + (0*2-5) + ... + (0*2-48) + (0*2-49) + (1*2-50) + (0*2-51) + (1*2-52)

Le résultat pour notre exemple est ainsi 0,6250000000000011102230246251565404236316680908203125

Combinaison exposant + mantisse

Une explication plus simple du fonctionnement de l’exposant combinée à la mantisse est de parler de fenêtres et de paliers. Imaginez que l’amplitude maximum permise par le type double précision soit découpée en fenêtres elles mêmes sous découpées en paliers avec :

  • 2045 fenêtres : 2047 exposants disponibles moins les exposants 0 et 2047 qui sont réservés.
    Une fenêtre est donc bornée entre deux puissances de 2 :
    • La première fenêtre (exposant 1) débute à 2-1022 (21-1023) et s’étend jusqu’à 2-1021.
    • La dernière fenêtre (exposant 2046) débute à 21023 (22046-1023) et s’étend jusqu’à 21024.
  • 252 paliers (4 503 599 627 370 496) par fenêtre : le nombre de valeurs de mantisse disponibles.
    Pour chaque fenêtre, un nombre de valeurs limité est ainsi disponible, la valeur d’incrémentation du palier (et donc la précision) dépend de la fenêtre à laquelle il se rattache.
    Le calcul de la précision est le suivant : (Taille de la fenêtre)/(nombre de valeurs de la mantisse)
    Par exemple (2-1021 - 2-1022) / 252 pour la première fenêtre.

Calcul de la mantisse et de l’exposant depuis une valeur réelle

Il est possible de calculer la mantisse et l’exposant d’une valeur réelle en suivant les étapes suivantes. la valeur d’exemple sera 1234,56789.

Étape 1 – Convertir le nombre en binaire :

1. Diviser la partie entière par 2 puis noter le reste (0 ou 1), répéter l’opération jusqu’à ce que l’opération soit égale à 0.
Le reste est noté en commençant par la droite.
Dans notre exemple, divisons successivement 1234 par 2 en récupérant le reste :

  1. 1234 ÷ 2 = 617 reste 0
  2. 617 ÷ 2 = 308 reste 1
  3. 308 ÷ 2 = 154 reste 0
  4. 154 ÷ 2 = 77 reste 0
  5. 77 ÷ 2 = 38 reste 1
  6. 38 ÷ 2 = 19 reste 0
  7. 19 ÷ 2 = 9 reste 1
  8. 9 ÷ 2 = 4 reste 1
  9. 4 ÷ 2 = 2 reste 0
  10. 2 ÷ 2 = 1 reste 0
  11. 1 ÷ 2 = 0 reste 1

Donc, 123410 ​= 100110100102

2. Multiplier la partie décimale par 2 puis noter la partie entière (0 ou 1), répéter l’opération un grand nombre de fois (voir l’étape 2 pour déterminer ce nombre).
La partie entière est notée en commençant par la gauche.
Dans notre exemple, multiplions successivement 0,56789 par 2, en prenant la partie entière de chaque produit :

  1. 0,56789 × 2 = 1,13578 (partie entière 1)
  2. 0,13578 × 2 = 0,27156 (partie entière 0)
  3. 0,27156 × 2 = 0,54312 (partie entière 0)
  4. 0,54312 × 2 = 1,08624 (partie entière 1)
  5. 0,08624 × 2 = 0,17248 (partie entière 0)
  6. 0,17248 × 2 = 0,34496 (partie entière 0)
  7. 0,34496 × 2 = 0,68992 (partie entière 0)
  8. 0,68992 × 2 = 1,37984 (partie entière 1)
  9. … on répète l’opération un nombre important de fois

En continuant, 0,5678910​ = 100100010110000100111101001100011011100111...2 (il s’agit d’une approximation)

La valeur totale approximée est ainsi : 10011010010,1001000101100001001111010011000110111001112

Étape 2 – Normaliser le nombre

Normaliser le nombre sous la forme 1,xyz x 2abc
Si le nombre débute par 0,x, il faudra alors utiliser une puissance négative.

La valeur normalisée doit avoir 52 décimales car la mantisse possède 52 bits. Il vous faudra ajuster le nombre de calculs nécessaires à l’étape 2. S’il ne vous est pas possible d’effectuer autant de calculs, complétez la mantisse par la valeur 0 (des bits nuls).

Pour notre exemple :
10011010010,1001000101100001001111010011000110111001112 = 1,0011010010100100010110000100111101001100011011100111 x 210

Dans le cas du nombre 0,0003, la valeur binaire est 0,000000000000001010001101101110001011101011000111000100001100101101.
La valeur normalisée est donc 1,010001101101110001011101011000111000100001100101101 x 2-14

Si votre valeur normalisée n’a pas 52 décimales, complétez avec des 0 : 1,0100110000000000000000000000000000000000000000000000 x 2y

Étape 3 – Extraire les informations

Les informations sont disponibles presque immédiatement :

  • L’exposant vous est donné par l’exposant de la valeur normalisée.
  • La mantisse est constituée de l’ensemble des valeurs suivant le « 1. ».

Dans notre exemple : 1,0011010010100100010110000100111101001100011011100111 x 210

  • L’exposant est 10
  • L’exposant décalé est 1033 (10 + 1023).
  • La fenêtre d’exposant est 210 à 211.
  • La valeur d’incrémentation et donc la précision est (211 - 210) / 252 = 2,27373675443232059478759765625 × 10-13.
  • La mantisse est : 0011010010100100010110000100111101001100011011100111
  • La valeur approchée stockée est 1234,567890000000033978722058236598968505859375 (calcul fait sur Wolfram Alpha).
    La valeur n’est précise qu’a 2,27 x 10-13 donc peut être arrondie à 1234,5678900000000 soit 1234,56789.

Valeurs spécifiques

La combinaison de valeurs spécifiques de la mantisse et de l’exposant permet de représenter certaines valeurs :

  • Exposant = 0 & mantisse = 0 : la valeur est 0
  • Exposant = 0 & mantisse != 0 : la valeur est dénormalisé
  • Exposant = 2047 & mantisse = 0 : la valeur est infini.
  • Exposant = 2047 & mantisse != 0 : NaN, la valeur n’est pas un nombre

Dénormalisation

La dénormalisation permet de stocker des nombres plus petits que 2-1022 (exposant le plus petit disponible). Pour cela, la méthode permettant de calculer la valeur finale change comme exposé plus haut :

  • Un nombre normalisé se calcule de la façon suivante :
    x = -1signe * (1 + mantisse) * 2exposant-1023
  • Alors qu’un nombre dénormalisé se calcule comme suivant :
    x = -1signe * (0 + mantisse) * 21-1023 = -1signe * 0,mantisse * 2-1022

Sans entrer dans les détails, le principe est de dire que dans ce mode de calcul, la valeur la plus haute de la mantisse (252) est alors égale à la valeur immédiatement avant la plus petite valeur en mode normalisé. Il devient ainsi possible d’augmenter l’amplitude.

D’un point de vue technique, le bit de décalage de la mantisse est alors défini à 0 et les programmes traitent la valeur de la mantisse comme dénormalisée.

Amplitude de valeurs

La combinaison de la mantisse et de l’exposant permet ainsi de calculer un ensemble de nombres avec une certaine précision. Cependant il n’est pas possible de représenter la totalité des nombres réels.

Une explication plus simple du fonctionnement est de parler de fenêtres et de paliers. La plage de valeurs utilisables est découpée en fenêtres elles mêmes sous découpées en paliers avec :

  • 2045 fenêtres : 2047 exposants disponibles moins les exposants 0 et 2047 réservés.
    Une fenêtre est donc bornée entre deux puissances de 2 :
    • La première fenêtre (exposant 1) débute à 2-1022 (21-1023) et s’étend jusqu’à 2-1021.
    • La dernière fenêtre (exposant 2046) débute à 21023 (22046-1023) et s’étend jusqu’à 21024.
  • 252 paliers (4 503 599 627 370 496) par fenêtre : le nombre de valeurs de mantisses disponibles.
    Pour chaque fenêtre, un nombre de valeurs limité est ainsi disponible, la valeur d’incrémentation du palier (et donc la précision) dépend de la fenêtre auquel le palier se rattache.
    Le calcul de la précision est le suivant : (Taille de la fenêtre)/(nombre de valeurs de la mantisse)
    Par exemple (2-1021 - 2-1022) / 252 pour la première fenêtre.

Partant de ce constat, voici un tableau de plusieurs valeurs intéressantes :

ExposantFenêtre d’exposantValeur d’incrémentation d’un palier
= Précision
MantisseValeur approchéeCommentaire
1de 2-1022 à 2-1021 : 4,940656458 x 10−32402,2250738585072013 x 10-308Plus petit nombre normalisé
1de 2-1022 à 2-10214,940656458 x 10−32412,2250738585072018 x 10−308Nombre immédiatement après le précédent
1023de 20 à 212,2204460492 × 10−1601,0Valeur 1
1023de 20 à 212,2204460492 × 10−1611,0000000000000002220446Nombre immédiatement après le précédent
1023de 20 à 212,2204460492 × 10−16450 359 962 737 0501,10000000000000008881784197001Plus proche palier de 1.1
1074de 251 à 2520,50
1075de 252 à 253104 503 599 627 370 496Le nombre suivant sera l’entier suivant
1076de 253 à 25420
2046de 21023 à 210241,99584030953 × 102921,79769313486231550856 x 10308Nombre immédiatement avant le plus grand nombre normalisé
2046de 21023 à 210241,99584030953 × 10292(1 + (1 − 2−52))1,79769313486231570814 x 10308Plus grand nombre normalisé

Ce tableau nous permet de voir plusieurs éléments :

  • tous les nombres ne peuvent être représentés exactement, c’est le cas de 1,1 qui est approximé à 1,100 000 000 000 000 088 817 841 970 01.
    Si un calcul est réalisé et qu’une valeur intermédiaire est trouvée, alors la valeur sera arrondie au palier le plus proche (voir la partie sur la précision).
  • Le passage d’une valeur à celle immédiatement après n’est pas régulier et dépend de l’ordre de grandeur d’un nombre.
  • Le point précédent implique que la précision d’un nombre dépend de son ordre de grandeur et peut aller de plusieurs centaines de chiffres après la virgule à plusieurs centaines de milliers en passant par un.

Les valeurs étant extrêmes, vous pouvez effectuer le calcul sur le site Wolfram Alpha.

Voici un tableau de quelques nombres dénormalisés :

ExposantMantisseValeur d’incrémentation d’un palier
= Précision
Valeur approchéeCommentaire
014,940656458 × 10-3244,940656458 × 10-324Plus petit nombre dénormalisé
024,940656458 × 10-3249,881 x 10-324Nombre dénormalisé suivant
02534,940656458 × 10-3242,2250738585072008 x 10-308Plus grand nombre dénormalisé

Ce tableau nous permet de voir que :

  • la valeur d’incrémentation des nombres dénormalisés est la même que le premier palier des nombres normalisé car le calcul se base sur le même exposant (la même fenêtre) : -1022 et la mantisse possède toujours le même nombre de valeurs (même nombre de paliers)
  • Le plus grand nombre dénormalisé est immédiatement avant la plus petite valeur normalisée ce qui permet une continuité des valeurs.

Voici enfin un dernier tableau présentant certains nombres spécifiques :

ExposantMantisseValeur approchéeCommentaire
000,0Zéro
20470InfiniInfini
20471NaNPremière valeur de NaN
2047253NaNDernière valeur de NaN

Précision et chiffres significatifs

La précision représente la plus petite différence pouvant être calculée entre deux valeurs consécutives. Cette précision varie selon l’ordre de grandeur de la valeur d’où le terme de « nombre à virgule flottante ».

Chaque valeur étant liée à une plage d’exposant (une fenêtre) et chaque plage d’exposant pouvant être découpée par 252 valeurs de mantisse (les paliers), la précision de chaque plage peut être estimée par l’approximation suivante :

précision (valeur / 252) (valeur / (4,5 x 1015))

Voici quelques ordres de grandeur de précision :

  • Plus petit nombre normalisé : 2,225 x 10-308 : 10-324
  • Valeur 1,0 : 10-16
  • Valeur 1,0 x 1016 : 1
  • Plus grand nombre normalisé : 1,798 x 10308 : 10292

Le terme de « nombre de chiffres significatifs » est également utilisé et représente le nombre de chiffre maximum qu’il faut pour percevoir une différence entre deux nombres en notation scientifique. Le nombre de chiffres significatifs est ici entre 15 et 16 car il faut au maximum 16 chiffres pour percevoir une différence entre deux nombres :

  • Plus petit nombre normalisé :
    • 2,225 073 858 507 201 3 x 10-308
    • 2,225 073 858 507 201 8 x 10−308
  • Valeur 1,0 :
    • 1,000 000 000 000 000 0 x 100
    • 1,000 000 000 000 000 2 x 100
  • Valeur dont l’incrémentation est de 1 à chaque palier :
    • 4,503 599 627 370 496 x 1015
    • 4,503 599 627 370 497 x 1016
  • Plus grand nombre normalisé : 1,798 x 10308 :
    • 1,797 693 134 862 315 5 x 10308
    • 1,797 693 134 862 315 7 x 10308

Dans le cas de PostgreSQL (et d’autres programmes) la valeur exacte stockée est arrondie selon la précision permise pour celle-ci (grâce à cette fonction dans le code source).
Ainsi, pour la valeur 1,1, le stockage en double précision donne la valeur exacte 1,10000000000000008881784197001 ce qui est précis à 10-29. Comme la précision permise pour l’ordre de grandeur de la valeur est de 2,22 x 10-16 près, la valeur renvoyée est arrondie selon cette précision donc 1,1000000000000000 soit 1,1.

Pour aller plus loin

Voici quelques ressources pour aller plus loin :

Précision sous PostGIS

Précision globale

Les coordonnées sous PostGIS sont stockées en utilisant le type double précision de PostgreSQL. La précision des valeurs dépend donc de l’ordre de grandeur et varie selon :

  • le système de projection.
  • la valeur de la coordonnée (et donc la position).
    En règle générale, il est possible d’appliquer la formule suivante pour avoir une idée de la précision possible : 10-16 x la valeur.

Voici quelques éléments de précision en lien avec différents systèmes de coordonnées :

WGS84 (EPSG:4326) :

  • Coordonnées X :
    • Amplitude : de -180 à 180
    • Précision : de 10-16 à 10-14 (sauf pour les coordonnées entre -1 et 1 qui peuvent atteindre une précision de 10-324).
  • Coordonnées Y :
    • Amplitude de -90 à 90
    • Précision : de 10-16 à 10-14 (sauf pour les coordonnées entre -1 et 1 qui peuvent atteindre une précision de 10-324).

WGS 84 / Pseudo-Mercator (EPSG:3857) :

  • Coordonnées X :
    • Amplitude : de -20037508,34 à 20037508,34
    • Précision : de 10-16 à 10-9 (sauf pour les coordonnées entre -1 et 1 qui peuvent atteindre une précision de 10-324).
  • Coordonnées Y :
    • Amplitude de -20048966,1 à 20048966,1
    • Précision : de 10-16 à 10-9 (sauf pour les coordonnées entre -1 et 1 qui peuvent atteindre une précision de 10-324).

Lambert 93 (EPSG:2154) :

  • Coordonnées X :
    • Amplitude : de -378305,81 à 1320649,57 (entre environ 100 000 et 1 100 000 pour la France métropolitaine)
    • Précision : de 10-16 à 10-10 (sauf pour les coordonnées entre -1 et 1 qui peuvent atteindre une précision de 10-324 mais qui ne sont presque jamais utilisées).
  • Coordonnées Y :
    • Amplitude de 6005281,2 à 7235612,72 (entre environ 6 033 000 et 7 130 000 pour la France métropolitaine)
    • Précision : 10-10

Pour le Lambert 93, la précision est entre 10-10 soit 0,1 nanomètres (1 ångström, la taille d’une molécule) et 10-16 soit 0,1 femtomètre (la taille d’un nucléon : proton ou neutron) ce qui suffit largement à positionner n’importe quoi…

Pour mieux visualiser la précision, il est possible de tester quelques requêtes dans PostGIS. Dirigeons nous vers Chamonix en Haute Savoie où se trouve un point possédant les coordonnées suivantes en Lambert 93 : X = 1000000 et Y = 6543210. Il s’agit de la terrasse d’un hôtel.

SELECT 
	id,
	geom_text,
	ST_AsText(geom_text::geometry(point, 2154)) AS geom
FROM (
	VALUES 
	(1, 'POINT (1000000 6543210)'),
	(2, 'POINT (1000000.0000000001 6543210.000000001)'),
	(3, 'POINT (1000000.00000000001 6543210.0000000001)'),
	(4, 'POINT (1000000.00000000008 6543210.0000000008)')
) AS t (id, geom_text)
;

Voici le résultat :

idgeom_textgeom
1POINT (1000000 6543210)POINT(1000000 6543210)
2POINT (1000000.0000000001 6543210.000000001)POINT(1000000.0000000001 6543210.000000001)
3POINT (1000000.00000000001 6543210.0000000001)POINT(1000000 6543210)
4POINT (1000000.00000000008 6543210.0000000008)POINT(1000000.0000000001 6543210.000000001)

Ce que l’on observe c’est que lorsqu’une coordonnée dont la précision est supérieure à ce qui est permis par le type double précision, alors celle-ci est automatiquement arrondie :

  • Dans le cas de la ligne 2, la valeur des coordonnées est permise par le type double précision, les coordonnées sont conservées telles quelles avec une précision « à la molécule près ».
  • Dans le cas de la ligne 3, les coordonnées X et Y possèdent une valeur dont l’ordre de grandeur est 106. La précision maximum est donc de 10-10 (106 x 10-16). Les valeurs sont donc arrondies à 10-10. Notez que pour la valeur de Y, la valeur de la partie entière impose un arrondi à 10-9 plutôt que 10-10.
  • Pour la ligne 4, le principe est le même mais l’arrondi se faisant au plus proche, nous récupérons la valeur au dessus.

Notez que la transformation des coordonnées dans le système WGS 84 (EPSG:4326) permet d’utiliser plus de décimales car les valeurs possèdent un ordre de grandeur plus faible :

SELECT 
	ST_AsText(
		ST_Transform(
			'POINT (1000000 6543210)'::geometry(point, 2154), 
			4326
		)
	) AS geom
;

Le résultat est POINT(6.87241320435217 45.92236395248839) : la précision est de 10-14 pour les coordonnées X et Y. Ici aussi vous pouvez voir que selon la valeur de la partie entière, la précision autorisée oscille entre 10-16 et 10-15 multipliée par l’ordre de grandeur.

Précision lors des calculs

Les calculs sous PostGIS sont réalisés en double précision et donc avec la précision détaillée précédemment. Par exemple, une distance entre deux objets aura une précision qui dépendra de l’ordre de grandeur de la valeur elle-même.

Lorsqu’il s’agit de comparer des géométries, le cas est plus complexe. Nous allons étudier la fonction ST_Intersects mais le principe est le même pour les autres fonctions de relation spatiale.

Relation sommet à sommet : il s’agit simplement de comparer les coordonnées entres-elles.

SELECT 
	id,
	geom_text,
	ST_AsText(geom_text::geometry(point, 2154)) AS geom,
	ST_Intersects(
		geom_text::geometry(point, 2154),
		'POINT(1000000 6543210)'::geometry(geometry, 2154)
	) AS st_intersects
FROM (
	VALUES 
	(1, 'POINT (1000000 6543210)'),
	(2, 'POINT (1000000.0000000001 6543210)'),
	(3, 'POINT (1000000.00000000001 6543210)'),
	(4, 'POINT (1000000.00000000008 6543210)'),
	(5, 'POINT (999999.99999999990 6543210)'),
	(6, 'POINT (999999.99999999997 6543210)')
) AS t (id, geom_text)
;
idgeom_textgeomst_intersects
1POINT (1000000 6543210)POINT(1000000 6543210)Oui
2POINT (1000000.0000000001 6543210.000000001)POINT(1000000.0000000001 6543210.000000001)Non
3POINT (1000000.00000000001 6543210.0000000001)POINT(1000000 6543210)Oui
4POINT (1000000.00000000008 6543210.0000000008)POINT(1000000.0000000001 6543210.000000001)Non
5POINT (999999.99999999990 6543210)POINT(999999.9999999999 6543210)Non
6POINT (999999.99999999997 6543210)POINT(1000000 6543210)Oui

Ainsi, la requête précédente permet de visualiser que la correspondance doit être exacte. Lorsqu’une coordonnées est « trop précise », elle est simplement arrondie afin de pouvoir être traitée puis comparée.

Notez que dans le cas du type geography, deux géométries s’intersectent si elles sont distantes de moins de 0,00001 mètre.

Relation sommet à segment : ce cas est plus complexe car il faut calculer la distance entre un sommet réel et un point virtuel situé sur le segment. Pour faire cela, plusieurs calculs sont réalisés, tous en double précision. Ainsi, la distance doit être inférieur à la précision permise par les coordonnées des objets comparés.

Ces opérations peuvent d’ailleurs entrainer des erreurs en raison de la nature imprécise des nombre en double précision. Il est donc possible d’avoir de faux négatifs : deux objets considérés comme ne s’intersectant pas alors qu’ils le sont ou inversement.

Voici une requête qui compare un point avec des lignes très proches de celui-ci :

SELECT 
	id,
	geom_line AS geom_text,
	ST_AsText(geom_line::geometry(linestring, 2154)),
	ST_Distance(
		geom_pt::geometry(point, 2154),
		geom_line::geometry(linestring, 2154)
	) AS st_distance,
	ST_Intersects(
		geom_pt::geometry(point, 2154),
		geom_line::geometry(linestring, 2154)
	) AS st_intersects
FROM (
	VALUES 
	(1, 'POINT (1000000 6543210)', 'LINESTRING(999999 6543210, 1000001 6543210)'),
	(2, 'POINT (1000000 6543210)', 'LINESTRING(999999 6543210, 1000001 6543210.0000000001)'),
	(3, 'POINT (1000000 6543210)', 'LINESTRING(999999 6543210, 1000001 6543210.000000001)'),
	(4, 'POINT (1000000 6543210)', 'LINESTRING(999999 6543210, 1000001 6543210.00000001)')
) AS t (id, geom_pt, geom_line)
;

Les objets sont alignés selon un axe globalement Est-Ouest. La variation de la coordonnée Y du sommet le plus à l’Est de la ligne permet de décaler légèrement celle-ci par rapport au point testé pour vérifier l’intersection des deux.

Le résultat est le suivant :

idgeom_textgeomst_distancest_intersects
1LINESTRING(999999 6543210, 1000001 6543210)LINESTRING (999999 6543210, 1000001 6543210)0.0Oui
2LINESTRING(999999 6543210, 1000001 6543210.0000000001)LINESTRING (999999 6543210, 1000001 6543210)0.0Oui
3LINESTRING(999999 6543210, 1000001 6543210.000000001)LINESTRING (999999 6543210, 1000001 6543210.000000001)0.0Non
4LINESTRING(999999 6543210, 1000001 6543210.00000001)LINESTRING (999999 6543210, 1000001 6543210.00000001)0.0000000055879354Non

La ligne 1 est une intersection parfaite. Vous pouvez observer que les géométries trop précises dans le cas de la ligne 2 sont ici aussi arrondies. En revanche, la ligne 3 met en valeur le fait que les arrondis de calcul entrainent une distance nulle alors que les objets ne s’intersectent pas ; en réalité, la distance est infinitésimale et est donc arrondie à 0. La ligne 4 montre un écart très petit.

Dernier exemple intéressant :

SELECT 
	id,
	geom_line AS geom_text,
	ST_AsText(geom_line::geometry(linestring, 2154)) AS geom_line,
	ST_Distance(
		geom_pt::geometry(point, 2154),
		geom_line::geometry(linestring, 2154)
	) AS st_distance,
	ST_Distance(
		geom_pt::geometry(point, 2154),
		geom_line::geometry(linestring, 2154)
	) * 1e100 AS st_distance_xe100,
	ST_Distance(
		geom_pt::geometry(point, 2154),
		geom_line::geometry(linestring, 2154)
	) * 1e308 AS st_distance_xe308,
	ST_Intersects(
		geom_pt::geometry(point, 2154),
		geom_line::geometry(linestring, 2154)
	) AS st_intersects,
	st_astext(
		ST_LineInterpolatePoint(
			geom_line::geometry(linestring, 2154),
			ST_LineLocatePoint(
				geom_line::geometry(linestring, 2154),
				geom_pt::geometry(point, 2154)
			)
		)
	) AS geom_point_proj,
	ST_Intersects(
		ST_LineInterpolatePoint(
			geom_line::geometry(linestring, 2154),
			ST_LineLocatePoint(
				geom_line::geometry(linestring, 2154),
				geom_pt::geometry(point, 2154)
			)
		),
		geom_line::geometry(linestring, 2154)
	) AS st_intersects_proj
FROM (
	VALUES 
	(1, 'POINT (1 0)', 'LINESTRING(0 0, 10 1e-100)'),
	(2, 'POINT (1 0)', 'LINESTRING(0 0, 10 1e-323)'),
	(3, 'POINT (1 0)', 'LINESTRING(0 0, 10 1e-324)')
) AS t (id, geom_pt, geom_line)
;

J’ai repris l’exemple d’une ligne avec un point mais l’ai amplifié : une ligne qui débute en X = 0 et finie en X = 10 est comparée à un point positionné en X = 2. Pour décaler la ligne du point, la fin de la ligne est déplacée selon l’axe Y d’une distance très petite (en limite du type double précision).

Je calcule la distance, ainsi que la distance multipliée par 10100 et 10300 pour éviter de passer sous la tolérance d’affichage des double précision de PostgreSQL.

Je détermine également un point projeté sur la ligne à partir du point testé.

idgeom_textgeom_linest_distancest_distance_xe100st_distance_xe300st_intersectsgeom_point_projst_intersects_proj
1LINESTRING(0 0, 10 0)LINESTRING(0 0, 10 0)000truePOINT(1 0)true
2LINESTRING(0 0, 10 1e-100)LINESTRING(0 0,10 1e-100)00.11099falsePOINT(1 1e-101)false
3LINESTRING(0 0, 10 1e-323)LINESTRING(0 0,10 1e-323)000falsePOINT(1 0)false
4LINESTRING(0 0, 10 1e-324)LINESTRING(0 0,10 0)000truePOINT(1 0)false

Vous pouvez voir que le point projeté n’intersecte la ligne sur laquelle il est censé être positionné que dans le premier cas. Dans le cas d’un décalage, même infinitésimale (dans le cas de la ligne 2, la distance entre la ligne et le point projeté est tellement infime qu’il n’est pas calculable avec PostGIS), il n’y a pas intersection.

Vous pouvez essayer avec des décalages plus importants (de l’ordre de 1 m par exemple) mais très souvent il n’y aura pas d’intersection car un point projeté sur un segment de ligne n’est qu’une approximation.

Pour avoir une intersection fiable, il doit y avoir une correspondance topologique : un sommet avec un sommet.

Réduire la précision

Comme avoir une cohérence topologique ne vous intéresse pas… vous voulez réduire la précision. Soit !

Il n’est pas possible de spécifier une précision à utiliser pour une base de données avec PostGIS. D’ailleurs le simple fait de vouloir contraindre la précision des coordonnées doit amener à réfléchir sur les propositions suivantes :

  • Il ne faut pas confondre précision de la mesure et précision des coordonnées : une mesure effectuée avec un GPS ayant une précision de ± 10 mètres doit être stockée telle que mesurée.
    Par exemple X = 33,5 m, alors la valeur réelle est comprise entre 23,5 et 43,5 m, réduire la précision à 10 mètres près entrainerait alors une perte d’information.
  • Lorsque les coordonnées d’un point sont calculées par interpolation, la valeur peut avoir plus de décimales que les valeurs utilisées dans l’interpolation.
    Par exemple, deux points distants dont les coordonnées sont au mètre près (pas de décimales) peuvent entrainer le calcul d’un point médian avec des coordonnées à décimales.
  • La réduction de la précision au niveau des coordonnées implique une perte de précision des valeurs issues de calculs à partir des coordonnées (aire, périmètre…).
  • La reprojection des coordonnées peut impliquer un changement du nombre de décimales.
    Par exemple, les coordonnées en WGS 84 (en degrés) impliquent un nombre plus important de décimales que les coordonnées en Lambert 93 (en mètres).

Il existe cependant plusieurs méthodes pour réduire la précision des coordonnées ou des calculs associés aux coordonnées.

Précision des coordonnées

ST_SnapToGrid : déplace les sommets de géométries en entrée sur une grille dont le pas est défini. Les points consécutifs déplacés au même endroit ainsi que les géométries nulles sont supprimés. Attention, des géométries invalides peuvent être créées.

ST_ReducePrecision : réduit la précision des coordonnées en les arrondissant selon la précision choisie. Les géométries créées sont toutes valides.

ST_QuantizeCoordinates : réduit la précision des coordonnées tout en opérant un traitement spécifique sur leur stockage afin d’en augmenter la compressibilité. Pour cela, les valeurs des bits non significatifs de la mantisse sont définies à zéro. Il s’agit de la meilleure méthode pour réduire le poids d’une géométrie.

Précision des calculs

Pour la précision des calculs, c’est différents. En effet, PostGIS utilise toujours la précision permise par le type double précision et seules quelques fonctions permettent d’utiliser une tolérance.

Voici les fonctions qui permettent de spécifier une tolérance (le paramètre gridsize) :

  • ST_Union
  • ST_UnaryUnion
  • ST_Difference
  • ST_Intersection
  • ST_SymDifference
  • ST_Subdivide

Cependant, dans la plupart des cas c’est la méthode qu’il faut adapter. Il s’agit ainsi d’utiliser une fonction réduisant la précision des coordonnées en conjonction d’une fonction classique ou d’utiliser une autre méthode que celle habituelle. Voici quelques cas et adaptations possibles.

ST_Intersects

Pour tester si deux géométries s’intersectent avec une certaine tolérance, il faut utiliser la fonction ST_DWithin() qui permet de déterminer si deux géométries sont situées à une certaine distance l’une de l’autre.

ST_Covers

Pour tester le recouvrement de deux géométries, vous pouvez appliquer une zone tampon à la géométrie « englobante » avec la fonction ST_Buffer(). Par exemple, une zone tampon de 5 mètres permet de savoir si une géométrie englobe une autre à 5 mètres près.

ST_Touches

Deux géométries se touchent si elle s’intersectent mais pas leurs « intérieurs ». Introduire une tolérance implique alors d’autoriser un léger chevauchement des intérieurs. Il faudra alors utiliser la fonction ST_DWithin pour valider que les géométries soient proches l’une de l’autre. Il faut ensuite valider que l’intersection des deux géométries soit minime en calculant l’aire d’intersection avec ST_Area(ST_Intersection()) et en la comparant à l’aire de chacune des géométries. Il conviendra de définir un pourcentage maximum (par exemple « moins de 1% d’intersection »).

Relation topologique : cas général

Pour le calcul des relations topologiques de manière générale, tout va dépendre du résultat voulu. Selon les cas, il faudra réduire la précision des coordonnées avant d’utiliser la fonction ou encore combiner plusieurs autres fonctions permettant de déterminer la relation voulue.

ST_Distance

Le calcul d’une distance n’est pas à proprement parlé lié à une tolérance mais la précision du résultat peut être modifié avec la fonction round().

Pour aller plus loin

Voici quelques liens pour approfondir le sujet :


Cet article vous a plu ?

N'hésitez pas à le partager, il interessera surement certains de vos contacts.

Les thèmes suivants contiennent des articles en lien avec celui-ci, allez faire un tour :

BDDPostGISPostgreSQLSIG double précisionexposantfloat8mantisseprécision

50%