[Postgis] Intersection raster vecteur

Version imprimable

Bonjour,

Dans une base Postgis, j'ai une couche raster et une couche vecteur (ligne).
Je voudrais récupérer les valeurs des pixels du raster qui intersectent ma ligne :
voici ma requête :
Code:

1 2 3 4 5 6 7 8 9 10 SELECT sampno, (ST_Intersection(r.no2, traj.geom)).geom road, (ST_Intersection(r.no2, traj.geom)).val elevation, r.rid from output_1_t100_27572 as r, (SELECT sampno, ST_Transform(geom, 27572) as geom from trajectories.corrected_gps_info_line_t where sampno='1108964') as traj where r.rid=27 and ST_Intersects(traj.geom, r.no2)
r est ma couche raster et traj ma couche vecteur.

Cette requête fonctionne, mais elle met plusieurs minutes à rendre le résultat !
Est-ce qu'il y a un moyen de l'écrire autrement pour optimiser le temps de réponse ? comme c'est pour afficher sur une interface cartographique web, il faudrait que le temps de réponse soit 'raisonnable'

Merci,
Nico

26/05/2017, 16h04
SQLpro

Avez-vous tenté de mettre un index spatial ?

Sinon il reste SQL Server bien plus rapide pour le traiteent des données spatial.
Pour info benchmark SIG SQL Server vs PostgreSQL: http://g-ernaelsten.developpez.com/t...-performances/

A +

Bonjour

J'essaie, mais quand je lance cette requete :

Code:

CREATE INDEX idx_output_1_t100_27572 ON output_1_t100_27572 USING GIST (no2)

J'ai le message d'erreur :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
 
ERREUR:  le type de données raster n'a pas de classe d'opérateurs par défaut pour la
méthode d'accès « gist »
HINT:  Vous devez spécifier une classe d'opérateur pour l'index ou définir une
classe d'opérateur par défaut pour le type de données.
********** Error **********
 
ERREUR: le type de données raster n'a pas de classe d'opérateurs par défaut pour la
méthode d'accès « gist »
SQL state: 42704
Hint: Vous devez spécifier une classe d'opérateur pour l'index ou définir une
classe d'opérateur par défaut pour le type de données.

Qu'est ce qu'il faut faire ??

Merci,
Nico

Salut
Je ne connais pas bien le SIG, mais qu'est ce qui empêche de remplacer la sous requête par la table correspondante?

Citation:

Envoyé par DiverSIG

Code:

1
2
3
4
5
6
7
8
9
10
SELECT sampno, 
(ST_Intersection(r.no2, traj.geom)).geom road, 
(ST_Intersection(r.no2, traj.geom)).val elevation,
r.rid
from output_1_t100_27572 as r,
	(SELECT sampno,  ST_Transform(geom, 27572) as geom
	from trajectories.corrected_gps_info_line_t
	where sampno='1108964') as traj
where r.rid=27
and ST_Intersects(traj.geom, r.no2)

Je propose donc...

Code:

1
2
3
4
5
6
7
SELECT sampno, 
(ST_Intersection(r.no2, ST_Transform(geom, 27572))).geom road, 
(ST_Intersection(r.no2, ST_Transform(geom, 27572))).val elevation,
r.rid
from output_1_t100_27572 as r,trajectories.corrected_gps_info_line_t as traj
where sampno='1108964' and r.rid=27
and ST_Intersects(traj.geom, r.no2)

29/05/2017, 09h55
DiverSIG

Bonjour,

Cette requête ne renvoie aucun résultat.
La requête avec la sous-requête renvoie le bon résultat (c'est (trop) long, mais ça marche)

Donc je suis plutôt dans la piste ''créer un index spatial'', mais je me débat avec mon message d'erreur (voir plus haut dans le post)

Merci,
Nico

Si j'essaie de lui spécifier un opérateur comme vu ici, j'ai le message d'erreur :

Code:

1
2
3
4
5
ERREUR:  la classe d'opérateur « gist_geometry_ops_2d » n'accepte pas le type de données raster
********** Error **********
 
ERREUR: la classe d'opérateur « gist_geometry_ops_2d » n'accepte pas le type de données raster
SQL state: 42804

Nico

30/05/2017, 01h35
alassanediakite

Salut

Citation:

Envoyé par DiverSIG

Cette requête ne renvoie aucun résultat...

Cela m'étonne énormément car je croie que les deux requêtes sont équivalentes sémantiquement. Peux-tu poster la requête que tu exécute et qui ne donne rien?
Pour l'index, une recherche sur google me donne...

Code:

CREATE INDEX nom_indexe ON nom_table USING GIST (ST_ConvexHull(colonne_raster))

@+
30/05/2017, 13h44
DiverSIG

J'ai exécuté texto la requête que tu as proposé : 0 réponse.
Pour la création d'index, j'ai déjà testé cela, mais j'ai une erreur de classe d'opérateur pour la méthode gist (cf. plus haut dans le post).

J'ai essayé de recharger ma couche avec les commandes raster2pgsql et psql, en spécifiant l'option -I dans le commande raster2pgsql (qui justement crée un index de type gist.
Ca fonctionne, dans pgadmin, je vois qu'il a crée un index nommé <nom de la couche>_st_convexhull_idx, donc c'est surement le même résultat que le create index que tu proposes.

Mais ça n'améliore pas le temps de réponse de la requête...

Nico

d'après ce que j'ai vu dans cet exemple, c'est pas gagné mon affaire !! (voir le chapitre 'intersecting the caribou buffers with the elevation rasters').

J'ai essayé de recharger ma couche en diminuant la taille des tuiles (20x20 au lieu de 100x100), j'ai maintenant une erreur de topologie :
Code:

1 2 3 4 5 6 7 8 9 ERREUR: GEOSIntersects: TopologyException: side location conflict at 255376.400374944 6252156.2600285504 CONTEXT: fonction PL/pgsql _st_intersects(geometry,raster,integer), ligne 22 à RETURN fonction PL/pgsql st_intersection(geometry,raster,integer), ligne 5 à affectation ********** Error ********** ERREUR: GEOSIntersects: TopologyException: side location conflict at 255376.400374944 6252156.2600285504 SQL state: XX000 Context: fonction PL/pgsql _st_intersects(geometry,raster,integer), ligne 22 à RETURN fonction PL/pgsql st_intersection(geometry,raster,integer), ligne 5 à affectation
Pourtant quand je fait une requete ST_IsValid sur ma sous-requête, il me dit geometrie OK.

Nico

Salut

Citation:

Envoyé par DiverSIG

J'ai exécuté texto la requête que tu as proposé : 0 réponse...

Tu as raison.
Ma requête doit être corriger au niveau du WHERE... traj.geom est à remplacer par ST_Transform(geom, 27572)
Code:

1 2 3 4 5 6 7 SELECT sampno, (ST_Intersection(r.no2, ST_Transform(geom, 27572))).geom road, (ST_Intersection(r.no2, ST_Transform(geom, 27572))).val elevation, r.rid from output_1_t100_27572 as r,trajectories.corrected_gps_info_line_t as traj where sampno='1108964' and r.rid=27 and ST_Intersects(ST_Transform(geom, 27572), r.no2)
Pour l'exemple cité, j'ai pas lu en entier, mais remarque la vétusté de produit...(PostgreSQL 8.4 et PostGIS 1.5.1)
Voici le lien que j'avais trouvé (copie/colle):http://mesange.educagri.fr/htdocs/si...x_spatial.html
@+

01/06/2017, 15h29
ced

Bonjour,

Une piste rapide (je n'ai pas mené de tests pour vérifier) : je pense que la non-utilisation des index, et donc la lenteur de la requête, vient en grande partie du fait que vous faite une intersection entre un raster et une couche vecteur reprojetée à la volée (utilisation de ST_Transform).
Du coup, à cause de cette reprojection, l'index GiST sur la couche vecteur n'est pas utilisable... Il faudrait nous afficher un EXPLAIN de la requête pour s'en assurer.

ced

1 pièce(s) jointe(s)

Bonjour,

J'ai donc transformé mon tracé GPS dans la même projection que mon raster (EPSG:3857).

La requête mouline toujours, si je l'arrête en cours d'exécution, j'ai le message :
Code:

1 2 3 4 5 6 7 8 ERREUR: Error performing intersection: InterruptedException: Interrupted! CONTEXT: fonction PL/pgsql st_intersection(geometry,raster,integer), ligne 10 à RETURN QUERY ********** Error ********** ERREUR: Error performing intersection: InterruptedException: Interrupted! SQL state: XX000 Context: fonction PL/pgsql st_intersection(geometry,raster,integer), ligne 10 à RETURN QUERY
si je lance un EXPLAIN de cette requête, j'ai le résultat :
Pièce jointe 281767

Il a l'air d'utiliser l'index sur le champ geom du raster...

Nico

02/06/2017, 12h54
ced
Deux remarques :
1. il faut mettre un index spatial sur la nouvelle couche vecteur (là, on voit un scan séquentiel) ;
2. il faut faire un ANALYZE corrected_gps_info_line_t_3857, parce que les statistiques de la table ne sont pas à jour.