jak nakładać plik shapefile i raster?

18

Mam plik kształtu z wielokątami. I mam globalny plik rastrowy. Chcę nałożyć wieloboki pliku kształtu na siatkę rastrową i obliczyć średnią wartość rastra dla każdego wielokąta.

Jak mogę to zrobić za pomocą GDAL, zapisując wyniki w pliku shapefile?

andreash
źródło
4
Czy GDAL jest jedynym narzędziem, którego chcesz użyć?
Simbamangu,
@Simbamangu nie, w zasadzie wszystko jest w porządku i byłoby wspaniale, gdyby było w Pythonie
andreash

Odpowiedzi:

9

W R możesz to zrobić

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

e jest wektorem ze średnią wartości komórek rastrowych dla każdego wielokąta.

RobertH
źródło
To nie jest R python, jak zadano w pytaniu
GM
6

Po otrzymaniu porady, którą znalazłem na liście mailingowej gdal-dev, użyłem StarSpan :

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

Wyniki są zapisywane w formacie CSV. W tym czasie było to już dla mnie wystarczające, ale powinno być możliwe wyodrębnienie pliku kształtu z tych informacji.

andreash
źródło
Wygląda na to, że StarSpan przeniósł się do GitHub. Zdobądź to tutaj .
Richard
5

Poniższy skrypt umożliwia wykonanie zadania za pomocą 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
źródło
4

Załaduj plik shapefile i raster do PostGIS 2.0 i wykonaj:

SELECT (ST_SummaryStats(ST_Clip(rast, geom))).*
FROM rastertable, geomtable
Pierre Racine
źródło
4

Nie sądzę, że GDAL jest najlepszym narzędziem do tego, ale możesz użyć gdal_rasterize, aby „wyczyścić” wszystkie wartości poza wielokątem.

Coś jak:

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

Program gdal_rasterize modyfikuje plik, dlatego tworzymy kopię do pracy. Zaznaczamy również określoną wartość (w tym przypadku zero), która ma być nodata. „-Burn 0 -b 1” oznacza wypalenie wartości zero w paśmie 1 pliku docelowego (work.tif). „-I” oznacza odwróconą rasteryzację, więc wypalamy wartości poza wielokątem zamiast wewnątrz niego. Komenda gdalinfo z -stats raportuje statystyki pasma. Wierzę, że wykluczy to wartość nodata (którą wcześniej oznaczyliśmy -a_nodata).

Frank Warmerdam
źródło
2

Przekształć plik kształtu w rastrze przez gdal_rasterize i użyj kodu w http://www.spatial-ecology.net/dokuwiki/doku.php?id=wiki:geo_tools do obliczenia statystyki strefowej dla każdego wielokąta. Możesz uruchomić http://km.fao.org/OFwiki/index.php/Oft-reclass, jeśli chcesz uzyskać tif ze statystykami rasters. Ciesz się kodem Ciao Giuseppe

Giuseppe
źródło
Czy zdarza ci się mieć kopię kodu, do którego się odwołujesz? Niestety link do pliku Python jest martwy.
ustroetz
1

Nie jest to możliwe przy użyciu GDAL. Możesz jednak użyć innych bezpłatnych narzędzi, np. Saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1
johanvdw
źródło
Poszedłem z tym podejściem, chociaż nazwa funkcji to tak naprawdę „Statystyka siatki dla wielokątów”.
bananowiec
1

Możesz także użyć rasterstats, który jest modułem Pythona zaprojektowanym do tego celu:

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

wprowadź opis zdjęcia tutaj

Następnie możesz uzyskać dostęp do atrybutu pierwszej strefy za pomocą:

mean_of_zone1 = listofzones[0]['mean']
GM
źródło
-2

możesz użyć narzędzia do obliczania statystyk punktów w arc gis, a narzędzie to można pobrać ze strony http://ianbroad.com/arcgis-toolbox-calculate-point-statistics-polygon-arcpy/

Suresh Goswami
źródło
2
„Narzędzie Oblicz statystyki punktów pobiera wejściową klasę wielokątów i punktów i wykorzystuje wybrane pole do znalezienia minimum, maksimum i średniej punktów i dodaje wyniki do elementu wielokąta”. ale to pytanie dotyczy klasy obiektów Wielokąt i Rastra, więc nie wydaje się odpowiednie.
PolyGeo