Difference between revisions of "Zonal Statistics or Extracting Raster Data"
(Set PropertyValue: TargetDate = 2015-01-20) |
|||
| (8 intermediate revisions by the same user not shown) | |||
| Line 24: | Line 24: | ||
| − | # Import a shapefile to R and use it to extract zonal statistics from a directory of rasters | + | <nowiki>#</nowiki> Import a shapefile to R and use it to extract zonal statistics from a directory of rasters |
| + | |||
| + | <nowiki>#</nowiki> Update working directory | ||
| − | |||
setwd("H:/Ian_GIS/Climate/precip/Tejon/Historic") | setwd("H:/Ian_GIS/Climate/precip/Tejon/Historic") | ||
| − | # Load raster package (may have to install first) | + | <nowiki>#</nowiki> Load raster package (may have to install first) |
| + | |||
library(raster) | library(raster) | ||
| − | # Make list of asci files in already specified working directory | + | <nowiki>#</nowiki> 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 | Clim_list = list.files(pattern = "mar.asc$") # "mar.asc$" queries files in directory ending in this pattern | ||
| − | # Check number of files in raster list | + | <nowiki>#</nowiki> Check number of files in raster list |
| + | |||
length = length(Clim_list) | length = length(Clim_list) | ||
| − | # Crete stack of rasters in list | + | <nowiki>#</nowiki> Crete stack of rasters in list |
| + | |||
Clim_stack = stack(Clim_list) | Clim_stack = stack(Clim_list) | ||
| − | # Create spatial dataset from shapefile (use points, polygons or lines) | + | <nowiki>#</nowiki> Create spatial dataset from shapefile (use points, polygons or lines) |
| − | library(maptools) | + | |
| + | library(maptools) | ||
| + | |||
points = readShapePoints("H:/Ian_GIS/Tejon_data/TejonPlots/TejonPoints") | points = readShapePoints("H:/Ian_GIS/Tejon_data/TejonPlots/TejonPoints") | ||
| + | |||
studyarea = readShapePoly("H:/Ian_GIS/Tejon_data/FullStudyArea/Tejon_study_area",IDvar="ID") | studyarea = readShapePoly("H:/Ian_GIS/Tejon_data/FullStudyArea/Tejon_study_area",IDvar="ID") | ||
| − | # Extract raster pixels from point dataset | + | <nowiki>#</nowiki> Extract raster pixels from point dataset |
for(j in 1:length) { | for(j in 1:length) { | ||
ex = extract(Clim_stack,points,method='simple',layer=1,nl=length, df=T) | ex = extract(Clim_stack,points,method='simple',layer=1,nl=length, df=T) | ||
ex$ID = NULL} | ex$ID = NULL} | ||
| − | # Extract study area mean from polygon dataset | + | <nowiki>#</nowiki> Extract study area mean from polygon dataset |
for(j in 1:length) { | for(j in 1:length) { | ||
bus = extract(Clim_stack,studyarea, fun=mean, layer=1, nl=length, df=T, na.rm=T) | bus = extract(Clim_stack,studyarea, fun=mean, layer=1, nl=length, df=T, na.rm=T) | ||
| Line 57: | Line 65: | ||
bus = t(bus) | bus = t(bus) | ||
| − | # Output csv file (point data) | + | <nowiki>#</nowiki> Output csv file (point data) |
| + | |||
write.csv(ex, "H:/Ian_GIS/Tejon/TejonPlots/Mar_precip_pts.csv", row.names=plots$Id, na="NA") | write.csv(ex, "H:/Ian_GIS/Tejon/TejonPlots/Mar_precip_pts.csv", row.names=plots$Id, na="NA") | ||
| − | # Output csv file (zonal polygon data) | + | <nowiki>#</nowiki> Output csv file (zonal polygon data) |
| + | |||
write.csv(bus, "H:/Ian_GIS/Tejon/FullStudyArea/Mar_precip_SA.csv", na="NA") | write.csv(bus, "H:/Ian_GIS/Tejon/FullStudyArea/Mar_precip_SA.csv", na="NA") | ||
| − | + | ||
<!-- Add any wiki Text above this Line --> | <!-- Add any wiki Text above this Line --> | ||
<!-- Do NOT Edit below this Line --> | <!-- Do NOT Edit below this Line --> | ||
{{#set: | {{#set: | ||
| − | |||
Expertise=R_programming| | Expertise=R_programming| | ||
Expertise=Raster| | Expertise=Raster| | ||
| + | Expertise=GIS| | ||
Owner=Ian_McCullough| | Owner=Ian_McCullough| | ||
Participants=Ian_McCullough| | Participants=Ian_McCullough| | ||
| + | Progress=100| | ||
| + | StartDate=2015-01-20| | ||
| + | TargetDate=2015-01-20| | ||
Type=Low}} | Type=Low}} | ||
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")
