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


	

Les GeoDataDays 2021 auront lieu les 15 et 16 septembre 2021 à Grenoble

 

Evénement national de référence et indépendant de la géographie numérique en France, cette rencontre est organisée par l'Afigéo et DécryptaGéo, en partenariat avec une plateforme régionale d'information géographique et des collectivités territoriales associées

Les inscriptions sont ouvertes et le programme disponible !

#1 Wed 05 May 2021 14:37

preliator
Membre
Date d'inscription: 17 Nov 2018
Messages: 382

Raster - statistique de zone avec PostGis

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)

Hors ligne

 

#2 Thu 06 May 2021 10:28

Nicolas Ribot
Moderateur
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1360

Re: Raster - statistique de zone avec PostGis

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

Hors ligne

 

#3 Thu 06 May 2021 10:46

preliator
Membre
Date d'inscription: 17 Nov 2018
Messages: 382

Re: Raster - statistique de zone avec PostGis

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 !

Hors ligne

 

#4 Thu 06 May 2021 10:56

Nicolas Ribot
Moderateur
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1360

Re: Raster - statistique de zone avec PostGis

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

Hors ligne

 

#5 Thu 06 May 2021 11:53

preliator
Membre
Date d'inscription: 17 Nov 2018
Messages: 382

Re: Raster - statistique de zone avec PostGis

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)

Hors ligne

 

Pied de page des forums

Powered by FluxBB

Partagez  |