bagaimana cara overlay shapefile dan raster?

17

Saya memiliki shapefile dengan poligon. Dan saya punya file raster global. Saya ingin overlay poligon shapefile ke grid raster dan menghitung nilai raster rata-rata untuk setiap poligon.

Bagaimana saya bisa melakukan ini menggunakan GDAL, menulis hasilnya ke shapefile?

andreash
sumber
4
Apakah GDAL satu-satunya alat yang ingin Anda gunakan?
Simbamangu
@Simbamangu tidak, pada dasarnya semuanya baik-baik saja, dan akan lebih bagus jika menggunakan Python
andreash

Jawaban:

9

Di R bisa Anda lakukan

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

e adalah vektor dengan rata-rata nilai sel raster untuk setiap poligon.

RobertH
sumber
Ini bukan python R seperti yang ditanyakan dalam pertanyaan
GM
6

Mengikuti saran yang saya dapat di milis gdal-dev, saya menggunakan StarSpan :

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

Hasilnya disimpan ke format CSV. Pada saat itu, itu sudah cukup bagi saya, tetapi harus dimungkinkan untuk memalsukan Shapefile dari info itu.

andreash
sumber
StarSpan tampaknya telah pindah ke GitHub. Dapatkan di sini .
Richard
4

Muat shapefile Anda dan raster Anda di PostGIS 2.0 dan lakukan:

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

Saya tidak berpikir GDAL adalah alat terbaik untuk ini, tetapi Anda dapat menggunakan gdal_rasterize untuk "menghapus" semua nilai di luar poligon.

Sesuatu seperti:

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 memodifikasi file, jadi kami membuat salinan untuk dikerjakan. Kami juga menandai beberapa nilai tertentu (nol dalam kasus ini) menjadi nodata. "-Burn 0 -b 1" berarti membakar nilai nol ke dalam band 1 dari file target (work.tif). "-I" berarti membalikkan rasterisasi sehingga kami membakar nilai di luar poligon bukan di dalamnya. Perintah gdalinfo dengan -stats melaporkan statistik band. Saya percaya ini akan mengecualikan nilai nodata (yang kami tandai sebelumnya dengan -a_nodata).

Frank Warmerdam
sumber
4

Script berikut memungkinkan Anda melakukan tugas dengan 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
sumber
2

Ubah file bentuk dalam raster oleh gdal_rasterize dan gunakan kode di http://www.spatial-ecology.net/dokuwiki/doku.php?id=wiki:geo_tools untuk menghitung statistik zona untuk setiap poligon. Anda dapat menjalankan http://km.fao.org/OFwiki/index.php/Oft-reclass jika Anda ingin mendapatkan tif dengan statistik raster Anda. Nikmati kode Ciao Giuseppe

Giuseppe
sumber
Apakah Anda memiliki salinan kode yang Anda rujuk? Sayangnya tautan ke file Python sudah mati.
ustroetz
1

Ini tidak mungkin menggunakan GDAL. Anda dapat menggunakan alat gratis lainnya, misalnya saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1
johanvdw
sumber
Saya menggunakan pendekatan ini, meskipun nama fungsinya sebenarnya adalah "Grid Statistics for Polygons".
bananafish
1

Anda juga dapat menggunakan rasterstats yang merupakan modul Python yang dirancang untuk tujuan ini:

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

masukkan deskripsi gambar di sini

Kemudian Anda dapat mengakses atribut zona pertama menggunakan:

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

Anda dapat menggunakan alat statistik penghitungan titik di arc gis dan alat ini dapat diunduh dari http://ianbroad.com/arcgis-toolbox-calculate-point-statistics-polygon-arcpy/

suresh Goswami
sumber
2
"Alat Calculate Point Statistics mengambil input Polygon dan kelas fitur Point dan menggunakan bidang yang dipilih untuk menemukan minimum, maksimum, dan rata-rata poin dan menambahkan hasilnya ke fitur poligon." tetapi pertanyaan ini adalah tentang kelas fitur Polygon dan Raster sehingga sepertinya tidak cocok.
PolyGeo