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

Suite à un problème technique intervenu entre le 22 et le 23 mars, nous avons du procéder dans la soirée du 25 mars, à la restauration de la base de données du 24 mars (matinée).

En clair, nous avons perdu vos contributions et inscriptions du dimanche 24 et du lundi 25 mars.
Nous vous prions de nous excuser.

#1 Wed 05 May 2021 14:37

preliator
Participant assidu
Date d'inscription: 17 Nov 2018
Messages: 433

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
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1534

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
Participant assidu
Date d'inscription: 17 Nov 2018
Messages: 433

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
Membre
Lieu: Toulouse
Date d'inscription: 9 Sep 2005
Messages: 1534

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
Participant assidu
Date d'inscription: 17 Nov 2018
Messages: 433

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