Le portail francophone de la géomatique

preliator · Wed 05 May 2021 14:37

Bonjour,

Je dispose d'une couche vecteur de polygones (40K) et d'un raster d'altitude que j'ai importé dans PostGis. Je souhaiterais obtenir la valeur moyenne d'altitude pour chaque polygone.

Tout d'abord, j'ai importé mon raster avec ma console :

Code:

raster2pgsql -s 2154 -C -I -M "K:\DATA\MNT\altitude.tif" -F -t 100x100 public.altitude_raster| psql -d db -U postgres -p 5432

Mes recherches m'ont menées vers la clause St_SummaryStats qui calcule un ensemble de statistiques (min, max, etc) dont ma chère moyenne.

Code:

SELECT  id, (stats).*
from (SELECT id, (ST_SummaryStats(a.rast, 1, true)) as stats
    FROM altitude_raster AS a, delete as b
    WHERE ST_Intersects(b.geom,a.rast)) as foo
order by 1

Bien que la requête fonctionne, je regarde des différences non négligeables avec l'outil "Statistique de zone" de QGis. Certaines fois, je suis à plusieurs centaines de mètres d'altitude de différence.

https://zupimages.net/viewer.php?id=21/18/wy24.jpg

J'ai également testé la requête suivante qui m'amène des résultats bien plus proches que sur QGis, mais avec l'apparition de valeurs nulles inexplicables.

Code:

select id, st_SummaryStats(ST_Union(ST_Clip(a.rast, 1, geom, true)))).*
from altitude_raster a, polyg b
where st_intersects(a.rast, b.geom)
group by 1

https://zupimages.net/viewer.php?id=21/18/jmb9.jpg

Bref, me suis-je trompé quelque part ?

Un grand merci.

Dernière modification par preliator (Thu 06 May 2021 10:45)

Nicolas Ribot · Thu 06 May 2021 10:28

Bonjour,

La condition st_intersects doit etre top généreuse: vous devez choper des tuiles qui touchent le pg mais ne l'intersecte pas, ou très peu.
Vous devriez faire l'intersection entre le raster et les pg et prendre en compte la part de chaque tuile contenue dans le polygone pour faire la moyenne.

Nicolas

preliator · Thu 06 May 2021 10:46

Hum, merci pour votre réponse.

Pour prendre en compte la part de chaque tuile, c'est le paramètre "1" dans

Code:

(ST_SummaryStats(a.rast, 1, true))

j'imagine ?

Merci !

Nicolas Ribot · Thu 06 May 2021 10:56

Non, plutot faire la moyenne des pixels qui intersectent, et surtout faire un clip entre le polygone et le raster.

Un bon exemple dans la doc Postgis: http://postgis.net/docs/RT_ST_SummaryStats.html
Dans la partie "Example: Summarize pixels that intersect buildings of interest"

La requete est plus complexe que la vôtre: d'abord l'intersection pg/raster, puis les stats par polygone.

Nicolas

preliator · Thu 06 May 2021 11:53

Merci je comprends mieux. Je me rappelle aussi que faire un St_clip sur mon raster me retourne des valeurs largement inférieures à 0. Je n'explique pas la raison. Voici le code lancé et le résultat.

Code:

(SELECT  id, (stats).*
    FROM (SELECT id, ST_SummaryStats(ST_Clip(rast, geom, true)) As stats
            FROM altitude_raster
            INNER JOIN pg
            ON ST_Intersects(pg.geom,rast)
        ) As foo
)

Et le résultat :

https://zupimages.net/viewer.php?id=21/18/iyhf.jpg

Dernière modification par preliator (Thu 06 May 2021 11:53)

Communauté	Ressources	L'emploi	Le marché	L'association	Fils RSS
Forum Blogs Wiki Planet Sites amis	Agenda Biblio Fiches AMO-CNIG Paris PCGIS Téléchargements	Offres d'emploi Banque de CV Formations Démarche Métiers	Geo-Entreprises Geo-Communiqué Appels d'offres	L'association Qui sommes-nous ? Nous aider Bienfaiteurs Partenaires Historique Statuts	Tout le site Job Géomatique Services GeoBlogs ...	Faire un don Mentions légales

Le portail francophone de la géomatique

#1 Wed 05 May 2021 14:37

Raster - statistique de zone avec PostGis

Code:

Code:

Code:

#2 Thu 06 May 2021 10:28

Re: Raster - statistique de zone avec PostGis

#3 Thu 06 May 2021 10:46

Re: Raster - statistique de zone avec PostGis

Code:

#4 Thu 06 May 2021 10:56

Re: Raster - statistique de zone avec PostGis

#5 Thu 06 May 2021 11:53

Re: Raster - statistique de zone avec PostGis

Code:

Pied de page des forums