¿Cómo superponer shapefile y raster?

17

Tengo un shapefile con polígonos. Y tengo un archivo ráster global. Quiero superponer los polígonos del archivo de forma en la cuadrícula de trama y calcular el valor medio de la trama para cada polígono.

¿Cómo puedo hacer esto usando GDAL, escribiendo los resultados en el shapefile?

andreash
fuente
44
¿Es GDAL la única herramienta que te gustaría usar?
Simbamangu
@Simbamangu no, básicamente todo está bien, y sería genial si estuviera en Python
andreash

Respuestas:

9

En R puedes hacer

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

e es un vector con la media de los valores de celda ráster para cada polígono.

RobertH
fuente
Esto es R, no Python como se hizo en la pregunta
GM
6

Siguiendo el consejo que recibí en la lista de correo de gdal-dev, utilicé StarSpan :

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

Los resultados se guardan en formato CSV. En ese momento, eso ya era suficiente para mí, pero debería ser posible de alguna manera falsificar un Shapefile a partir de esa información.

andreash
fuente
StarSpan parece haberse mudado a GitHub. Consíguelo aquí .
Richard
4

Cargue su archivo de forma y su ráster en PostGIS 2.0 y haga:

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

No creo que GDAL sea la mejor herramienta para esto, pero puede usar gdal_rasterize para "borrar" todos los valores fuera del polígono.

Algo como:

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

El programa gdal_rasterize modifica el archivo, por lo que hacemos una copia para trabajar. También marcamos algún valor particular (cero en este caso) como nodata. El "-burn 0 -b 1" significa grabar un valor de cero en la banda 1 del archivo de destino (work.tif). El "-i" significa rasterización invertida, por lo que quemamos valores fuera del polígono en lugar de dentro de él. El comando gdalinfo con -stats informa sobre estadísticas de banda. Creo que excluirá el valor de nodata (que marcamos anteriormente con -a_nodata).

Frank Warmerdam
fuente
4

El siguiente script le permite realizar la tarea con 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
fuente
2

Transforme el archivo de forma en ráster mediante gdal_rasterize y use el código en http://www.spatial-ecology.net/dokuwiki/doku.php?id=wiki:geo_tools para calcular la estadística zonal para cada polígono. Puede ejecutar http://km.fao.org/OFwiki/index.php/Oft-reclass si desea obtener un tif con la estadística de sus rásteres. Disfruta el código Ciao Giuseppe

Giuseppe
fuente
¿Tienes una copia del código al que te refieres? Lamentablemente, el enlace al archivo Python está muerto.
ustroetz
1

Esto no es posible con GDAL. Sin embargo, podría usar otras herramientas gratuitas, por ejemplo, saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1
johanvdw
fuente
Seguí este enfoque, aunque el nombre de la función es en realidad "Estadísticas de cuadrícula para polígonos".
bananafish
1

También puede usar rasterstats que es un módulo de Python diseñado para este propósito:

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

ingrese la descripción de la imagen aquí

Luego puede acceder al atributo de la primera zona usando:

mean_of_zone1 = listofzones[0]['mean']
GM
fuente
-2

puede usar la herramienta de cálculo de estadísticas de puntos en arc gis y esta herramienta se puede descargar desde http://ianbroad.com/arcgis-toolbox-calculate-point-statistics-polygon-arcpy/

suresh Goswami
fuente
2
"La herramienta Calcular estadísticas de puntos toma una clase de entidad de Polígono y Punto de entrada y usa un campo seleccionado para encontrar el mínimo, máximo y promedio de los puntos y agrega los resultados a la entidad de polígono". pero esta pregunta se trata de una clase de entidad Polígono y un Ráster, por lo que parece poco probable que sea adecuada.
PolyGeo