Le portail francophone de la géomatique

Bonjour,

Vous pouvez montrer la requete SQL, la creation des index et le plan de la requete ?
Analyse a bien ete lancé apres la creation des index ?

Nicolas

Fais juste un explain si la requete est longue/ne finit pas (explain analyze va lancer la requête en plus du plan)

Nico

Merci,
Pour le plan, ca serait mieux la sortie texte de pg, plus facile a lire je trouve.

Concernant la creation de l'index spatial, je n'ai jamais utilisé les options que tu indiques.
Tu as essayé une creation simple, style: 

CREATE INDEX ON public.elevationp USING gist(geom);
analyze elevationp;

Pour voir si ca change le plan ?

Nicolas

Ce sont des geom 3D ?
Meme si c'est le cas, tu peux indexer en 2D je pense

Sinon ce qui marche très vite pour récuperer des valeurs dans un raster à partir de vecteurs, c'est d'utiliser gdallocationinfo sur le raster direct: ca peut choper des centaines de milliers de valeurs par seconde et les stocker dans postgis.

J'ai pas mal joué avec ca ces temps-ci et ca permet de faire des croisements massifs entre vecteurs et raster (style des linéaires sur des pays entiers vs des rasters a 30m de réso).

Le principe est :
• de créer une grille virtuelle correspondant au raster (un pg de la grille = 1 pixel du raster) UNIQUEMENT sous les geometries à traiter (les batiments par ex). Ca va vite en général, sauf pour des super gros volumes, style batiments de france vs raster a 1m par ex.
Pour les ~ 110 000 batiments Paris sur une grille 30m, ca prend qq secondes de faire cette grille.

• de prendre les centroids de cette grille, de les passer a gdallocationinfo qui renvoie les valeurs dans pg direct (ultra rapide). C'est faisable en une passe avec un peu de shell, style:

psql -c "select x, y from gridtable order by row_id" | gdallocationinfo -valonly -geoloc dem.tif | psql -c "COPY gridtableval(val) from STDIN"

gridtableval est la table de résultat à créer avant qui permet de stocker la sortie de gdallocationinfo.
En ajoutant une pk generated always as identity, on peut lier cette table de résultat sur la table gridtable pour avoir la valeur raster de chaque pixel de la grille

• de faire le SQL qui donne la synthèse par batiment (ici le max d'alti pour pixels de chaque batiment).

• S'il faut faire des analyses plus poussées, par ex prendre en compte la partie de pixel qui intersecte la geom, on travaille sur ces pixels uniquement et pas tous les pixels du raster.

Avec cette méthode, des process qui prenaient des heures prennent qq minutes maintenant.

Nico

Bonjour,

Et merci à tous les deux,
Je vais faire différents tests un peu plus tard ...
La méthode de Nicolas me plaît bien, car c'est effectivement la question que je creuse, utilisation des fonctionnalités raster(in db) de postgis.
Comme l'IGN a libéré plein de données raster ça peut permettre des analyses pointues.

Ici calculer l'élévation maximale thérique d'un bâtiment.
Ce qui m'attire aussi dans les fonctionnalités raster c'est la possibilité de créer mes propres raster avec des valeurs dont je maîtrise la valeur. Particulièrement dans le cadre de calculs ou fusion entre rasters de différentes résolutions.

J'ajoute que combiner les résultats avec l'OrthoHR avec des RASTEROP ça ouvre des pistes prometteuses.

Je vous tiens au courant de la suite.

"3.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1"

Re,

Pour l'index en fait il ne faut pas préciser d'opérateur :

CREATE INDEX IF NOT EXISTS gix_b
    ON cadastre.batiment USING gist
    (the_geom)
    WITH (buffering=auto)
    TABLESPACE pg_default;

Là l'index est pris en compte par le query planner

Hello,

Cool que ca marche.

La solution que j'évoquais avec gdallocationinfo ne se base pas sur postgis raster, mais juste postgis pour les vecteurs, et gdal pour accéder aux rasters.
Le trick étant alors coté postgis de créer rapidement une grille "virtuelle" représentant le raster pour avoir les coord des pixels a choper avec gdal.
Dans mon cas, c'était pour éviter le chargement de gros raster dans la BD alors que peu d'info de pixels étaient nécessaires par rapport au raster entier.

La méthode postgis_raster marche bien aussi comme tu dis Christophe: là on a la totale: les valeurs rasters, chaque pixel en tant que polygone, et des index spatiaux pour aller vite. On peut en plus fabriquer n'importe quel nouveau raster par requête, ce qui est très puissant.

Dans ton croisement vecteur/raster, tu peux aller plus direct je pense: la CTE p1 ne sert pas de ce que je comprends.

Nicolas