Difference between revisions of "Zonal Statistics or Extracting Raster Data"

Latest revision as of 04:23, 21 January 2015

Zonal Statistics

Let's say you have a series of annual precipitation grids for a given site. Within this site you have a set of study plots and you want to extract precipitation values from each grid for each of your study plots. This is great application of zonal statistics. If you plots are points, zonal statistics will simply extract precipitation values for the location of each point. If your plots are lines or polygons, zonal statistics can be used to calculate the mean (or some other statistic) precipitation for each plot.

Skill Relevance and Usefulness

Use of zonal statistics is widely applicable across disciplines that utilize spatial datasets.

What You'll Need

The approach shown here uses the mpatools and raster libraries in R http://cran.r-project.org/web/packages/maptools/maptools.pdf http://cran.r-project.org/web/packages/raster/raster.pdf

A surprising number of people is unaware that you can use GIS shapefiles in R.The maptools library imports .shp points, lines and polygons

Example Code

This is a code that you could adapt for your own purposes; hopefully it is enough to guide you. You'll obviously need to designate your own directories and name your own files. Also, this code contains lines for polygons and points; you may very well want one or the other. The list function might be confusing. I used pattern-matching to pick out only March precipitation grids from my directory (which contained grids for all months). The pattern="mar.asc$" tells R to pick out files with a specific type of name, in this case ascii grids for March.

This example extracts March precipitation data for a set of study plots in the Tehachapi Mountains, CA. The outputs are csv files written to a user-specified directory.

# Import a shapefile to R and use it to extract zonal statistics from a directory of rasters

# Update working directory

setwd("H:/Ian_GIS/Climate/precip/Tejon/Historic")

# Load raster package (may have to install first)

library(raster)

# Make list of asci files in already specified working directory

Clim_list = list.files(pattern = "mar.asc$") # "mar.asc$" queries files in directory ending in this pattern

# Check number of files in raster list

length = length(Clim_list)

# Crete stack of rasters in list

Clim_stack = stack(Clim_list)

# Create spatial dataset from shapefile (use points, polygons or lines)

library(maptools)

points = readShapePoints("H:/Ian_GIS/Tejon_data/TejonPlots/TejonPoints")

studyarea = readShapePoly("H:/Ian_GIS/Tejon_data/FullStudyArea/Tejon_study_area",IDvar="ID")

# Extract raster pixels from point dataset for(j in 1:length) {

 ex = extract(Clim_stack,points,method='simple',layer=1,nl=length, df=T)
 ex$ID = NULL}

# Extract study area mean from polygon dataset for(j in 1:length) {

 bus = extract(Clim_stack,studyarea, fun=mean, layer=1, nl=length, df=T, na.rm=T)
 bus$ID = NULL}

bus = t(bus)

# Output csv file (point data)

write.csv(ex, "H:/Ian_GIS/Tejon/TejonPlots/Mar_precip_pts.csv", row.names=plots$Id, na="NA")

# Output csv file (zonal polygon data)

write.csv(bus, "H:/Ian_GIS/Tejon/FullStudyArea/Mar_precip_SA.csv", na="NA")