how to overlay shapefile and raster?

Question

I have a shapefile with polygons. And I have a global raster file. I want to overlay the shapefile's polygons onto the raster grid and calculate the mean raster value for each polygon.

How can I do this using GDAL, writing the results to the shapefile?

Answer 1

In R you can do

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

e is a vector with the mean of the raster cell values for each polygon.

Answer 2

Load your shapefile and your raster in PostGIS 2.0 and do:

SELECT (ST_SummaryStats(ST_Clip(rast, geom))).*
FROM rastertable, geomtable

Answer 3

I don't think GDAL is the best tool for this, but you can use gdal_rasterize to "clear" all the values outside the polygon.

Something like:

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

The gdal_rasterize program modifies the file, so we make a copy to work on. We also mark some particular value (zero in this case) to be nodata. The "-burn 0 -b 1" means burn a value of zero into band 1 of the target file (work.tif). The "-i" means invert rasterization so we burn values outside the polygon instead of inside it. The gdalinfo command with -stats reports on band statistics. I believe it will exclude the nodata value (which we marked earlier with -a_nodata).

Answer 4

Following advice I got on the gdal-dev mailing list, I used StarSpan:

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

The results are saved to CSV format. At that time, that was already enough for me, but it should be possible to somehow forge a Shapefile from that info.

Answer 5

Transform the shape file in raster by gdal_rasterize and the use the code in http://www.spatial-ecology.net/dokuwiki/doku.php?id=wiki:geo_tools to calculate zonal statistic for each polygon. You can run http://km.fao.org/OFwiki/index.php/Oft-reclass if you want obtain a tif with your rasters statistic. Enjoy the code Ciao Giuseppe

Answer 6

This is not possible using GDAL. You could use other free tools however, eg saga gis:

saga_cmd shapes_grid "Grid Values to Shapes" -GRIDS=grid.sgrd -POLYGONS=in.shp -SHAPES=out.shp-NODATA -TYPE=1