Discussion :

[DiverSIG]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[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 +

[DiverSIG]

Bonjour

J'essaie, mais quand je lance cette requete :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

CREATE INDEX idx_output_1_t100_27572 ON output_1_t100_27572 USING GIST (no2)

J'ai le message d'erreur :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[alassanediakite]

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?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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)

@+

[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

[DiverSIG]

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

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[alassanediakite]

Salut

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

@+

[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

[DiverSIG]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[alassanediakite]

Salut

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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
@+

[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

[DiverSIG]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :


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

Nico

[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.