Few days back a friend asked me to help him with a task. He wanted to clip a raster with a polygon vector file and then calculate mean of clipped raster pixel values.

I figured out that it could be done in following steps

Load raster and vector

Read vector and get all polygons

While Looping though the polygons clip raster based on the selected polygon calculate mean write in a csv with polygon’s unique key



I used Tropical Rainfall Measuring Mission (TRMM, January 1998) raster image and African boundaries shape file for this script (see below).

I decided to write a script in R. I looked at some of the examples on stagexchange.com and came across an answer (link)which helped in writing the code below (almost 80% of the code)

# install.packages("raster") #into your R console library(raster) #Note: # (1) I assume that raster and vector are in one EPSG + vector is overlaping raster (polygons) # (2) Vector and raster are in in project directory #read shape file shp = shapefile("Africa_boundaries.shp") # load reaster img = raster("rasterfile.tif"); # filter by unique polygons ( this will go into loop) country <- c("Algeria","Angola","Bahrain","Benin","Botswana","Burkina Faso","Burundi","Cameroon","Cape Verde","Central African Republic","Chad","Comoros","Congo","Djibouti","Egypt","Equatorial Guinea","Eritrea","Ethiopia","Gabon","Gambia, The","Gaza Strip","Ghana","Guinea","Guinea-Bissau","Ivory Coast","Kenya","Lebanon","Lesotho","Liberia","Libya","Madagascar","Malawi","Mali","Mauritania","Mauritius","Morocco","Mozambique","Namibia","Niger","Nigeria","Reunion","Rwanda","Sao Tome and Principe","Senegal","Seychelles","Sierra Leone","Somalia","South Africa","Sudan","Swaziland","Tanzania, United Republic of","Togo","Tunisia","Uganda","Western Sahara","Zaire","Zambia","Zimbabwe"); output <- matrix(ncol=3, nrow=7) for (i in 1:length(country)) { # now filter shape file to get one polygon record filtered = shp[match(toupper(country[i]),toupper(shp$name)),] #now use the mask function to clip based on polygon record rr <- mask(img, filtered) # get mean of all cells country_sum <-cellStats(rr, mean, na.rm=TRUE); # write record in csv output[i,1] <- i; output[i,2] <- country[i]; output[i,3] <- country_sum; } write.csv(output, file ="data.csv")

I should have looped through shape file to read each polygon rather than using a variable like country. Next time I will do that :).

The extraction and cellStats function took little bit of time. Finally, I joined my csv with shape file (using qGIS, see below).

Its better to know more than one programming language and if you are looking for an other language R is a strong candidate especially for GIS/RS.

reference:

http://gis.stackexchange.com/questions/61243/clipping-a-raster-in-r