Nous utilisons des cookies pour vous garantir la meilleure expérience sur notre site. Si vous continuez à utiliser ce dernier, nous considèrerons que vous acceptez l'utilisation des cookies. J'ai compris ! ou En savoir plus !.
banniere

Le portail francophone de la géomatique


Toujours pas inscrit ? Mot de passe oublié ?
Nom d'utilisateur    Mot de passe              Toujours pas inscrit ?   Mot de passe oublié ?

Annonce

Pour sa 21ème année, l’association GeoRezo a toujours besoin de vous !

10€ = 1 mois de frais bancaires ; 15€ = 12 mois de nom de domaine ; 30€ = 1 semaine de location des serveurs …

Faire un don 

Retrouver nos membres bienfaiteurs

#1 Thu 05 October 2023 08:51

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Soucis d'index spatial

Bonjour,

Je possède une table de polygones (bâtiments), et une table de polygones issus d'un raster (4000*6000) soit 24 millions de polygones vs une centaines dans la table bâtiments.

Après création d'un idex gist sur les deux tables quand je tente une requête d'intersection avec un WHERE st_intersects(geoma,geomb) le query planneur (lol) refuse d'utiliser l'index spatial ...

I Wonder Why ?

Bonne journée.


Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#2 Thu 05 October 2023 09:13

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

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

Hors ligne

 

#3 Thu 05 October 2023 10:18

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Bonjour,

Merci Nicolas.
Non j'ai pas fait analyse après la création. je teste.


Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#4 Thu 05 October 2023 10:20

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

ChristopheV a écrit:

Bonjour,

Merci Nicolas.
Non j'ai pas fait analyse après la création. je teste.


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

Nico

Hors ligne

 

#5 Thu 05 October 2023 10:24

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Re

Analyse lancé sur la table, pas de changement

Création index table issue de raster

-- Index: gist_p10

-- DROP INDEX IF EXISTS public.gist_p10;

CREATE INDEX IF NOT EXISTS gist_p10
    ON public.elevationp USING gist
    (geom gist_geometry_ops_nd)
    WITH (buffering=auto)
    TABLESPACE pg_default;

PJ Query Explain

Query
SELECT idbatiment,dur,max(elevation),the_geom
FROM elevationp,cadastre.batiment
WHERE st_intersects(geom,the_geom) --AND idbatiment>=2789 AND idbatiment<=2790
GROUP BY idbatiment,dur,the_geom

Dernière modification par ChristopheV (Thu 05 October 2023 10:25)


Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#6 Thu 05 October 2023 10:47

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

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: 

Code:

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

Pour voir si ca change le plan ?

Nicolas

Hors ligne

 

#7 Thu 05 October 2023 10:51

tumasgiu
Membre
Lieu: Ajaccio
Date d'inscription: 5 Jul 2010
Messages: 1161

Re: Soucis d'index spatial

Hello,

est ce que tu as fait un vacuum analyze sur les deux tables ?

Hors ligne

 

#8 Thu 05 October 2023 10:57

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

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

Hors ligne

 

#9 Thu 05 October 2023 10:59

tumasgiu
Membre
Lieu: Ajaccio
Date d'inscription: 5 Jul 2010
Messages: 1161

Re: Soucis d'index spatial

Je me souviens aussi que dans de vieilles versions de postgis, st_intersects ne provoquait pas forcement l'utilisation de l'index et qu'il fallait rajouter la condition d'intersection de bounding box a.geom && b.geom

Hors ligne

 

#10 Thu 05 October 2023 11:15

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

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:

Code:

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

Hors ligne

 

#11 Thu 05 October 2023 11:17

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

tumasgiu a écrit:

Je me souviens aussi que dans de vieilles versions de postgis, st_intersects ne provoquait pas forcement l'utilisation de l'index et qu'il fallait rajouter la condition d'intersection de bounding box a.geom && b.geom


Oui au fait, que donne postgis_full_version() ?


(je vous recommande pg16 + postgis 3.4 (geos 12): ca marche bien !)

Hors ligne

 

#12 Thu 05 October 2023 11:30

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Bonjour,

Et merci à tous les deux,
Je vais faire différents tests un peu plus tard ... smile
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"


Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#13 Sat 07 October 2023 07:06

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Bonjour,

La piste de Nicolas est la bonne.

Un extrait de code qui permet d'obtenir très rapidement la "vectorisation" d'un raster en carrés et points, plus valeur du raster.
Ici le raster est nommé "elevation". Nota il est découpé en tuile de 100*100, conservation de l'index de la tuile (rid).

Code:

with p as (SELECT rid,rast,st_squaregrid(1,st_envelope(rast)) rg FROM elevation
LIMIT 100000)

SELECT rid,st_value(rast,st_centroid((rg).geom)) val_rast,st_astext((rg).geom)as sgeom,
st_astext(st_centroid((rg).geom)) as pgeom
FROM p
ORDER BY rid DESC

100 000 pixels tratés en 2s ...

Merci Nicolas.


Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#14 Sat 07 October 2023 07:49

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Re,

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

Code:

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


Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#15 Sat 07 October 2023 09:56

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Bonjour,

La solution finale.

Code:

with p as (SELECT idbatiment,dur,rid,the_geom ,rast
FROM elevation,cadastre.batiment
WHERE st_intersects(the_geom,st_envelope(rast)) AND NOT st_bandisnodata(rast)),

p1 as (SELECT idbatiment,dur,rid,rast,the_geom as rid_bat_geom
FROM p 
ORDER BY rid),
p2 as (SELECT rid,st_centroid((st_squaregrid(1,st_envelope(rast))).geom) pg FROM elevation)

SELECT idbatiment,dur,p2.rid,st_value(rast,pg),pg 
--INTO elevation_batiment
FROM p1,p2 
WHERE st_intersects(rid_bat_geom,pg) AND p1.rid=p2.rid

Ce qui donne un temps d'exécution de 37.35 secondes pour un raster de 4 000 * 6 000
et 4 161 bâtiments dans l'emprise du raster.

Et un rendu sympathique cf. PJ

NB : Ceci permet aussi de faire du pur SQL sans passer par gdgdallocationinfo()


Fichier(s) joint(s) :
Pour accéder aux fichiers vous devez vous inscrire.

Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

#16 Tue 10 October 2023 11:48

Nicolas Ribot
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1554

Re: Soucis d'index spatial

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

Hors ligne

 

#17 Tue 10 October 2023 17:19

ChristopheV
Membre
Lieu: Ajaccio
Date d'inscription: 7 Sep 2005
Messages: 3199
Site web

Re: Soucis d'index spatial

Salut,

Oui Nicolas j'ai pas bien compris, vire la CTE et regarde la différence du query planner, du coup j'ai laissé me disant que cela allait forcément interpeller quelques uns ...

Je me permets le raster in DB car j'ai une BD par projet (BBox). Donc quelques gros raster gênent pas ...

L'idée étant qu'à partir d'une BD centrale regroupant l'ensemble des données d'une région (ou département), je créé un script qui génère l'ensemble des données dont j'ai besoin dans la base du projet, plus quelques pages de textes (de .teX en fait) avec graphiques et tableaux.Les MNT de l'IGN par exemple, je les laisse en fichiers jpg2000 avec du out-db pour avoir les BBox et chemin et je génère uniquement pour les besoins du projet en in db, et/ou avec un st_asgdalraster() + lo_export() pour génèrer directement les jpeg2000 du projet ...

Comme cela pour chaque études j'ai une BD spatiale prête et un corps d'analyse pré-écrit.

Cela permet de gagner du temps et de se consacrer à la tâche essentielle l'analyse, quand la vision carto et statistique est prête automatiquement ça facilite les choses.

Nota : L'algo précédent peut-être amélioré en tuilant le raster plus fortement, comme la comparaison ce fait au niveau des tuiles, et qu'il y a un index sur le rid ...


Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close

Hors ligne

 

Pied de page des forums

Powered by FluxBB