Wie überlagere ich Shapefile und Raster?

17

Ich habe ein Shapefile mit Polygonen. Und ich habe eine globale Rasterdatei. Ich möchte die Polygone des Shapefiles auf das Rasterraster legen und den mittleren Rasterwert für jedes Polygon berechnen.

Wie kann ich das mit GDAL machen und die Ergebnisse in das Shapefile schreiben?

andreash
quelle
4
Ist GDAL das einzige Tool, das Sie verwenden möchten?
Simbamangu
@Simbamangu nein, im Grunde ist alles in Ordnung, und es wäre toll, wenn es in Python wäre
andreash

Antworten:

9

In R können Sie tun

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

e ist ein Vektor mit dem Mittelwert der Rasterzellenwerte für jedes Polygon.

RobertH
quelle
Dies ist R nicht Python wie in der Frage gestellt
GM
6

Nach Beratung ich auf der gdal-dev Mailing - Liste bekam, habe ich StarSpan :

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

Die Ergebnisse werden im CSV-Format gespeichert. Zu diesem Zeitpunkt war das schon genug für mich, aber es sollte irgendwie möglich sein, ein Shapefile aus diesen Informationen zu fälschen.

andreash
quelle
StarSpan scheint zu GitHub gewechselt zu sein. Hol es dir hier .
Richard
4

Laden Sie Ihr Shapefile und Ihr Raster in PostGIS 2.0 und führen Sie Folgendes aus:

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

Ich denke nicht, dass GDAL das beste Werkzeug dafür ist, aber Sie können gdal_rasterize verwenden, um alle Werte außerhalb des Polygons zu "löschen".

Etwas wie:

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

Das Programm gdal_rasterize ändert die Datei, sodass wir eine Kopie erstellen, an der gearbeitet werden kann. Wir markieren auch einen bestimmten Wert (in diesem Fall Null) als Nodata. "-Burn 0 -b 1" bedeutet, dass ein Wert von Null in Band 1 der Zieldatei (work.tif) gebrannt wird. Das "-i" bedeutet invertierte Rasterisierung, sodass Werte außerhalb des Polygons anstatt innerhalb desselben gespeichert werden. Der Befehl gdalinfo mit -stats gibt Auskunft über die Bandstatistik. Ich glaube, es wird den Nodata-Wert ausschließen (den wir zuvor mit -a_nodata markiert haben).

Frank Warmerdam
quelle
4

Mit dem folgenden Skript können Sie die Aufgabe mit GDAL ausführen: 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
quelle
2

Transformieren Sie die Formdatei in Raster mit gdal_rasterize und verwenden Sie den Code in http://www.spatial-ecology.net/dokuwiki/doku.php?id=wiki:geo_tools , um die Zonenstatistik für jedes Polygon zu berechnen. Sie können http://km.fao.org/OFwiki/index.php/Oft-reclass ausführen, wenn Sie ein TIF mit Ihrer Raster-Statistik erhalten möchten. Genießen Sie den Code Ciao Giuseppe

Giuseppe
quelle
Haben Sie zufällig eine Kopie des Codes, auf den Sie sich beziehen? Leider ist der Link zur Python-Datei nicht mehr vorhanden.
Uströtz
1

Dies ist mit GDAL nicht möglich. Sie können jedoch auch andere kostenlose Tools verwenden, z. B. saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1
johanvdw
quelle
Ich habe diesen Ansatz gewählt, obwohl der Funktionsname eigentlich "Grid Statistics for Polygons" ist.
Bananafish
1

Sie können auch Rasterstats verwenden, da es sich um ein Python-Modul handelt, das für diesen Zweck entwickelt wurde:

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

Bildbeschreibung hier eingeben

Dann können Sie auf das Attribut der ersten Zone zugreifen, indem Sie Folgendes verwenden:

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

Sie können das Werkzeug zur Berechnung der Punktstatistik in arc gis verwenden. Dieses Werkzeug kann unter http://ianbroad.com/arcgis-toolbox-calculate-point-statistics-polygon-arcpy/ heruntergeladen werden.

suresh Goswami
quelle
2
"Das Werkzeug Punktstatistik berechnen verwendet eine Eingabe-Polygon- und Punkt-Feature-Class und verwendet ein ausgewähltes Feld, um das Minimum, Maximum und den Durchschnitt der Punkte zu ermitteln und die Ergebnisse zum Polygon-Feature hinzuzufügen." Bei dieser Frage handelt es sich jedoch um eine Polygon-Feature-Class und ein Raster, sodass es unwahrscheinlich ist, dass sie geeignet ist.
PolyGeo