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 …
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
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)
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
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 ...
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()
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