comment superposer un fichier de formes et un raster?

18

J'ai un fichier de formes avec des polygones. Et j'ai un fichier raster global. Je veux superposer les polygones du fichier de formes sur la grille raster et calculer la valeur raster moyenne pour chaque polygone.

Comment puis-je faire cela en utilisant GDAL, en écrivant les résultats dans le fichier de formes?

andreash
la source
4
GDAL est-il le seul outil que vous souhaitez utiliser?
Simbamangu
@Simbamangu non, fondamentalement tout va bien, et ce serait génial si c'était en Python
andreash

Réponses:

9

Dans R, vous pouvez faire

library(raster)
library(rgdal)
r <- raster('raster_filename')
p <- readOGR('shp_path', 'shp_file')
e <- extract(r, p, fun=mean)

e est un vecteur avec la moyenne des valeurs des cellules raster pour chaque polygone.

RobertH
la source
Ce n'est pas R python comme demandé dans la question
GM
6

Suite aux conseils que j'ai reçus sur la liste de diffusion gdal-dev, j'ai utilisé StarSpan :

starspan --vector V --raster R1 R2 ... --stats mystats.csv avg mode

Les résultats sont enregistrés au format CSV. À ce moment-là, cela me suffisait déjà, mais il devrait être possible de créer un Shapefile à partir de ces informations.

andreash
la source
StarSpan semble avoir migré vers GitHub. Obtenez-le ici .
Richard
5

Le script suivant vous permet d'effectuer la tâche avec GDAL: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#calculate-zonal-statistics

# Calculates statistics (mean) on values of a raster within the zones of an polygon shapefile

import gdal, ogr, osr, numpy

def zonal_stats(input_value_raster, input_zone_polygon):

    # Open data
    raster = gdal.Open(input_value_raster)
    driver = ogr.GetDriverByName('ESRI Shapefile')
    shp = driver.Open(input_zone_polygon)
    lyr = shp.GetLayer()

    # get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # reproject geometry to same projection as raster
    sourceSR = lyr.GetSpatialRef()
    targetSR = osr.SpatialReference()
    targetSR.ImportFromWkt(raster.GetProjectionRef())
    coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()
    geom.Transform(coordTrans)

    # Get extent of geometry
    ring = geom.GetGeometryRef(0)
    numpoints = ring.GetPointCount()
    pointsX = []; pointsY = []
    for p in range(numpoints):
            lon, lat, z = ring.GetPoint(p)
            pointsX.append(lon)
            pointsY.append(lat)
    xmin = min(pointsX)
    xmax = max(pointsX)
    ymin = min(pointsY)
    ymax = max(pointsY)

    # Specify offset and rows and columns to read
    xoff = int((xmin - xOrigin)/pixelWidth)
    yoff = int((yOrigin - ymax)/pixelWidth)
    xcount = int((xmax - xmin)/pixelWidth)+1
    ycount = int((ymax - ymin)/pixelWidth)+1

    # create memory target raster
    target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, gdal.GDT_Byte)
    target_ds.SetGeoTransform((
        xmin, pixelWidth, 0,
        ymax, 0, pixelHeight,
    ))

    # create for target raster the same projection as for the value raster
    raster_srs = osr.SpatialReference()
    raster_srs.ImportFromWkt(raster.GetProjectionRef())
    target_ds.SetProjection(raster_srs.ExportToWkt())

    # rasterize zone polygon to raster
    gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])

    # read raster as arrays
    banddataraster = raster.GetRasterBand(1)
    dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float)

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float)

    # mask zone of raster
    zoneraster = numpy.ma.masked_array(dataraster,  numpy.logical_not(datamask))

    # calculate mean of zonal raster
    return numpy.mean(zoneraster)
ustroetz
la source
4

Chargez votre fichier de formes et votre raster dans PostGIS 2.0 et faites:

SELECT (ST_SummaryStats(ST_Clip(rast, geom))).*
FROM rastertable, geomtable
Pierre Racine
la source
4

Je ne pense pas que GDAL soit le meilleur outil pour cela, mais vous pouvez utiliser gdal_rasterize pour "effacer" toutes les valeurs en dehors du polygone.

Quelque chose comme:

gdal_translate -a_nodata 0 original.tif work.tif
gdal_rasterize -burn 0 -b 1 -i work.tif yourpolygon.shp -l yourpolygon
gdalinfo -stats work.tif
rm work.tif

Le programme gdal_rasterize modifie le fichier, nous faisons donc une copie sur laquelle travailler. Nous marquons également une valeur particulière (zéro dans ce cas) comme nodata. Le "-burn 0 -b 1" signifie graver une valeur de zéro dans la bande 1 du fichier cible (work.tif). Le "-i" signifie une tramage inversée, donc nous brûlons les valeurs à l' extérieur du polygone au lieu de l'intérieur. La commande gdalinfo avec -stats rend compte des statistiques de bande. Je pense que cela exclura la valeur nodata (que nous avons marquée plus tôt avec -a_nodata).

Frank Warmerdam
la source
2

Transformez le fichier de forme en raster par gdal_rasterize et utilisez le code dans http://www.spatial-ecology.net/dokuwiki/doku.php?id=wiki:geo_tools pour calculer les statistiques zonales pour chaque polygone. Vous pouvez exécuter http://km.fao.org/OFwiki/index.php/Oft-reclass si vous souhaitez obtenir un tif avec vos statistiques de rasters. Profitez du code Ciao Giuseppe

Giuseppe
la source
Avez-vous une copie du code auquel vous faites référence? Malheureusement, le lien vers le fichier Python est mort.
ustroetz
1

Ce n'est pas possible avec GDAL. Vous pouvez cependant utiliser d'autres outils gratuits, par exemple les saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1
johanvdw
la source
Je suis allé avec cette approche, bien que le nom de la fonction soit en fait "Grid Statistics for Polygons".
bananafish
1

Vous pouvez également utiliser rasterstats qui est un module Python conçu à cet effet:

from rasterstats import zonal_stats
listofzones = zonal_stats("polygons.shp", "elevation.tif",
            stats="mean")

entrez la description de l'image ici

Ensuite, vous pouvez accéder à l'attribut de la première zone en utilisant:

mean_of_zone1 = listofzones[0]['mean']
GM
la source
-2

vous pouvez utiliser l'outil de calcul de statistiques de points dans arc gis et cet outil peut être téléchargé à partir de http://ianbroad.com/arcgis-toolbox-calculate-point-statistics-polygon-arcpy/

suresh Goswami
la source
2
"L'outil Calculer les statistiques de points prend une classe d'entités polygone et point en entrée et utilise un champ sélectionné pour trouver le minimum, le maximum et la moyenne des points et ajoute les résultats à l'entité surfacique." mais cette question concerne une classe d'entités polygones et un raster, il semble donc peu probable qu'elle convienne.
PolyGeo