¿Cómo superponer shapefile y raster?

14

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

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

    
pregunta andreash 18.03.2012 - 13:32

8 respuestas

6

Siguiendo consejos entré en la lista de correo de gdal-dev , Usé 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 forjar un Shapefile a partir de esa información.

    
respondido por el andreash 20.06.2012 - 18:30
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 celdas ráster para cada polígono.

    
respondido por el RobertH 22.03.2012 - 21:05
4

Cargue su shapefile y su raster en PostGIS 2.0 y haga:

SELECT (ST_SummaryStats(ST_Clip(rast, geom))).*
FROM rastertable, geomtable
    
respondido por el Pierre Racine 18.03.2012 - 21:11
4

No creo que GDAL sea la mejor herramienta para esto, pero puedes 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) para ser nodata. "-Burn 0 -b 1" significa grabar un valor de cero en la banda 1 del archivo de destino (work.tif). La "-i" significa rasterización inversa, por lo que quemamos los valores fuera del polígono en lugar de dentro de él. El comando gdalinfo con -stats informa sobre las estadísticas de banda. Creo que excluirá el valor de nodata (que marcamos anteriormente con -a_nodata).

    
respondido por el Frank Warmerdam 17.04.2012 - 00:13
4

El siguiente script le permite realizar la tarea con GDAL: enlace

# 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)
    
respondido por el ustroetz 29.05.2013 - 01:19
2

Transforme el archivo de formas en ráster por gdal_rasterize y use el código en enlace para calcular la estadística zonal para cada polígono. Puede ejecutar enlace si desea obtener un tif con su estadística de rasters. Disfruta el codigo Ciao Giuseppe

    
respondido por el Giuseppe 23.08.2012 - 00:41
1

Esto no es posible usando GDAL. Puede 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
    
respondido por el johanvdw 18.03.2012 - 20:05
-2

puede utilizar la herramienta para calcular estadísticas de puntos en arc gis y esta herramienta se puede descargar desde enlace

    
respondido por el suresh Goswami 14.10.2016 - 06:19

Lea otras preguntas en las etiquetas