R Spatial

Report
Spatial
1
Types of Spatial Data
R Spatial packages
S4 R objects
Projections
Spatial methods
•
•
•
Outline
SpatialPoints
SpatialLines/SpatialPolygons
SpatialPixels
SpatialPoints - Sample Dataset
•
•
•
Importing/Exporting
Spatial information and slots
Displaying
SpatialPolygons - Sample Dataset
•
•
•
•
Importing/Exporting
Spatial information and slots
Displaying
Reprojecting
SpatialPixels - Sample Dataset
•
•
•
•
•
•
Importing/Exporting rasters
Spatial information and slots
Displaying
Deriving
Thematic rasters (importing, spatial info, reclassing)
Temporary files
2
Spatial Data in R
## Types of Spatial Data
Points - a set of single point locations
Lines
- an ordered set of points, connected by straight line segments
Polygons – an area, with enclosed lines, possibly containing holes
Rasters – a collection of points or rectangular cells, in grid format
## R Spatial packages
sp
Basic R classes for handling geospatial data
maptools Read and write shapefiles. Cannot read the projection file.
rgeos
Interface to spatial geometry for sp objects
rgdal
Supports GDAL raster formats and OGR vector formats..
Retains projection information when reading and writing (PROJ.4 library)
raster
Spatial data analysis for rasters.
## Help links for spatial tools in R
# Summary of tools for Spatial data analysis
http://cran.r-project.org/web/views/Spatial.html
3
Spatial Data in R
# Set working directory
path <- "W:/Techniques/Projects/Peru/Rworkshop"
setwd(path)
# Load libraries
library(sp)
library(maptools)
library(rgeos)
library(rgdal)
library(raster)
4
Spatial Data in R
## Spatial classes (sp package)
getClass("Spatial")
Bivand, et al. 20085
Spatial Data in R
## S4 R objects
New-style classes with formal definitions that specify the name and type of the
object's components, called slots.
## The Spatial class has 2 pre-defined slots
1. bbox – The bounding box, consisting of a matrix of coordinates defining the
extent of the spatial object.
2. proj4string – The class object defining the coordinate reference system (CRS)
CRS contains datum and projection
datum – a surface that represents the shape of the earth
projection – renders the surface of an ellipsoid as a plane
6
Spatial Data in R
## Projections
# rgdal package
PROJ.4 - Cartographic Projections Library
A library of projection functions to perform respective forward and inverse
transformation of cartographic data.
# Format of projection string (proj4string)
"+proj=prj +ellps=GRS80 +datum=datum +no_defs”
# List of common projections
http://www.remotesensing.org/geotiff/proj_list/
# More proj4string options..
" +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"
" +proj=utm +zone=13 +datum=NAD83"
" +proj=utm +ellps=WGS84 +datum=WGS84 +zone=11, +units=m +towgs84=0,0,0"
" +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m"
" +no_defs +ellps=GRS80 +towgs84=0,0,0"
7
Spatial Data in R
## Common attributes of R objects:
# method
- a function associated with a particular type of object (ex. summary, print)
## Spatial methods (sp package)
Bivand, et al. 20088
SpatialPoints
Bivand, et al. 20089
SpatialLines/Polygons
Lists of List Objects
List of Objects
Object
Bivand, et al. 200810
SpatialPixels
Similar to SpatialPoints with regularly spaced coordinates.
Bivand, et al. 200811
Sample Dataset
Utah, USA
12
Sample Dataset
Uinta Mountains,
Utah,USA
Highest East-West oriented
mountain range in the
contiguous U.S. - up to
13,528 ft (4,123 m)
High Uinta Wilderness
## Vegetation
5 different life zones:
1. shrub-montane
2. aspen
3. lodgepole pine
4. spruce-fir
5. alpine
13
SpatialPoints
14
Sample Dataset
SpatialPoints - Importing/Exporting
## Import data frame with X/Y coordinates
##
(FIA plot data with fuzzed/swapped coordinates)
plt <- read.csv("PlotData/plt.csv", header=TRUE,
stringsAsFactors=FALSE)
## Generate SpatialPoints
x <- "LON"
y <- "LAT"
prj <- "longlat"
datum <- "NAD83"
prj4str <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"
ptshp <- SpatialPointsDataFrame(plt[,c(x,y)], plt,
proj4string = CRS(prj4str))
class(ptshp)
# Export SpatialPoints
help(writeOGR)
writeOGR(ptshp, "Outfolder", "ptshp", driver="ESRI Shapefile")
15
Sample Dataset
SpatialPoints – Spatial Info
## Get spatial information for ptshp
summary(ptshp)
mode(ptshp)
class(ptshp)
# General info and summary statistics for attribute
# data of ptshp
# Get mode of ptshp
# Get class of ptshp
is.projected(ptshp)
proj4string(ptshp)
# Check if ptshp has projection information
# Set or get projection information
str(ptshp)
# Structure of ptshp
16
Sample Dataset
SpatialPoints - Slots
## Get slot names and access slots for SpatialPoints
slotNames(ptshp)
[email protected]
[email protected]
[email protected]
[email protected]
#
#
#
#
#
Get
Get
Get
Get
Get
name of ptshp components (slots)
data frame attributes of ptshp
coordinates of ptshp
extent of ptshp
projection information
head([email protected])
dim(ptshp)
# Get first 6 rows of data frame of ptshp
# Get dimensions of data frame of ptshp
17
Sample Dataset
SpatialPoints - Displaying
## Display ptshp
plot(ptshp)
# Display ptshp points with red dots, half size
plot(ptshp, col="RED", pch=16, cex=0.5)
# Display subsets of ptshp
plot(ptshp[ptshp$INVYR==2002,], col="BLUE")
plot(subset(ptshp, INVYR==2003), col="RED", add=TRUE)
# pch codes
# Colors
colors()
## Help for point graphics
http://127.0.0.1:11789/library/graphics/html/points.html
18
Sample Dataset
Exercise
## Exercise 1
## 1.1 Get the min and max coordinates of the boundary of ptshp.
## 1.2 Display points where elevation (ELEVM) are >= 3000 in dark green.
## 1.3 Add points where elevation (ELEVM) are < 3000 in yellow.
## 1.4 Overlay the point with the maximum elevation (ELEVM) in red.
19
SpatialPolygons
20
Sample Dataset
SpatialPolygons
# SpatialPolygons – file names
# Area of Interest (AOI) boundary (Uinta Mountains, UT)
aoifn <- "uintaN_aoi.shp"
# National Forest System boundary within AOI
nfsfn <- "uintaN_nfs.shp"
# High Uinta Wilderness boundary within AOI
wildfn <- "uintaN_wild.shp"
21
Sample Dataset
SpatialPolygons - Importing/Exporting
## Importing shapefiles using OGR drivers (rgdal)
## OGR functions
From the Geospatial Data Abstraction Library (GDAL), interfaced through rgdal.
Reads spatial data, including the spatial reference.
readOGR – reads OGR data source and layer into an R Spatial object
writeOGR – writes data out using supported drivers
Note: There are 2 main arguments.. that may take different forms, depending on the
available drivers.
dsn – data source name
layer – the name of the layer
Note: A driver is a software component loaded on demand to enable a device to work with
a computer's operating system.
ogrDrivers()
help(readOGR)
22
Sample Dataset
SpatialPolygons - Import
## Importing polygon shapefiles
## (AOI boundary, NFS boundary, Wilderness boundary)
# dsn – data source name (folder with layers)
# layer – the name of the layer (no extension)
## Set dsn
dsn <- "SpatialData"
## Spatial layer names
aoinm <- "uintaN_aoi"
nfsnm <- "uintaN_nfs"
wildnm <- "uintaN_wild"
## Import shapefiles
bndpoly <- readOGR(dsn=dsn, layer=aoinm, stringsAsFactors=FALSE)
nfspoly <- readOGR(dsn=dsn, layer=nfsnm, stringsAsFactors=FALSE)
wildpoly <- readOGR(dsn=dsn, layer=wildnm, stringsAsFactors=FALSE)
23
Sample Dataset
SpatialPolygons – Spatial Info
## Get spatial information for bndpoly
bndpoly
summary(bndpoly)
mode(bndpoly)
typeof(bndpoly)
class(bndpoly)
proj4string(bndpoly)
bbox(bndpoly)
dim(bndpoly)
#
#
#
#
#
#
#
#
Get general information of polygon
Summary information for polygon
Get mode of bndpoly
Get typeof
Get class of bndpoly
Set or get projection information of bndpoly
Get extent of bndpoly
Get dimension of bndpoly
str(bndpoly)
# Get structure of bndpoly
24
Sample Dataset
SpatialPolygons - slots
(1)
(2)
(3)
# Get slot names for class SpatialPolygons (1)
getClass("SpatialPolygons")
slotNames(bndpoly)
@polygons
@plotOrder
@bbox
@proj4string
@data
# the list of Polygons objects
# the order to plot the Polygons
# bounding box slot
# coordinate reference system slot
# associated data frame (for class SpatialPolygonsDataFrame objects)
25
Sample Dataset
SpatialPolygons – Displaying
## Display bndpoly
plot(bndpoly)
# Display polygon shapefile
plot(bndpoly, col="blue")
# Display polygon with blue fill color
plot(bndpoly, border="red", lwd=3) # red outline and line width=3
plot(bndpoly, col="blue", border="red", lwd=3) # blue fill
plot(bndpoly, col="blue", border="red", lwd=3, bg="black", axes=TRUE)
# with background colored and axes labels
26
Sample Dataset
SpatialPolygons – Spatial Info
# Get spatial information for nfspoly
nfspoly
summary(nfspoly)
dim(nfspoly)
# Get general information of nfspoly
# Summary information for nfspoly
# Dimensions of nfspoly
str(nfspoly)
# Get structure of nfspoly
27
Sample Dataset
SpatialPolygons - slots
(1)
(2)
(3)
# Get slot names for class SpatialPolygons (1)
getClass("SpatialPolygons")
slotNames(nfspoly)
@polygons
@plotOrder
@bbox
@proj4string
@data
# the list of Polygons objects
# the order to plot the Polygons
# bounding box slot
# coordinate reference system slot
# associated data frame (for class SpatialPolygonsDataFrame objects)
28
Sample Dataset
SpatialPolygons - slots
(1)
(2)
(3)
# Get slot names for class Polygons (2)
getClass("Polygons")
slotNames([email protected][[1]])
@polygons[[1]]@Polygons
@polygons[[1]]@plotOrder
@polygons[[1]]@labpt
@polygons[[1]]@ID
@polygons[[1]]@area
# First feature
# a list of rings (islands/holes) that make up the feature
# the order to plot the feature
# label point coordinates of the feature
# unique identifier of the feature
# area of the feature (in units of feature projection)
29
Sample Dataset
SpatialPolygons - slots
(1)
(2)
(3)
# Get slot names for class Polygon (3)
getClass("Polygon")
slotNames(slot([email protected][[1]], "Polygons")[[1]])
@polygons[[1]]@Polygons[[1]]
@polygons[[1]]@Polygons[[1]]@labpt
@polygons[[1]]@Polygons[[1]]@area
@polygons[[1]]@Polygons[[1]]@coords
# first ring of feature - class 'Polygon'
# label point coordinates of the feature ring
# area of the feature ring
# coordinates of the feature ring
30
Sample Dataset
SpatialPolygons - slots
# Accessing slots of SpatialPolygons
[email protected]
[email protected]
[email protected]
[email protected]
#
#
#
#
Data frame attributes of shp
Get info of each polygon within Spatial Polygons
Get order of polygons within Spatial Polygons
Get extent of bndpoly
# Accessing slots of Polygons
slot(nfspoly, "polygons")
slotsPolygons <- slot(nfspoly, "polygons")
typeof(slotsPolygons)
length(slotsPolygons)
# (2)slots
# (2)slots
sapply(slotsPolygons, function(x)slot(x, "ID"))
sapply(slotsPolygons, function(x)slot(x, "area"))
# (2)slots
# (2)slots
slotsPolygon <- slot([email protected][[1]], "Polygons")
slotsPolygon <- slot([email protected][[2]], "Polygons")
# (3)slots
# (3)slots
31
Sample Dataset
SpatialPolygons - Displaying
## Display nfspoly and wildpoly
plot(bndpoly, col="dark blue")
plot(nfspoly, border="yellow", add=TRUE, lwd=2)
plot(wildpoly, add=TRUE, border="green")
# Add points – in red
plot(ptshp, add=TRUE, col="red", pch=16, cex=0.5)
32
Sample Dataset
SpatialPolygons - Reprojecting
## Notice, the point layer did not display. Check the projections.
projection(bndpoly)
projection(ptshp)
## Reproject ptshp to projection of bndpoly
polyprj <- projection(bndpoly)
ptshpprj <- spTransform(ptshp, CRSobj=CRS(polyprj))
projection(ptshpprj)
# Add points – in red
plot(ptshpprj, add=TRUE, col="red", pch=16, cex=0.5)
# Add thicker boundary line
plot(bndpoly, lwd=3, add=TRUE)
# Add labels for nfspoly
text(nfspoly, nfspoly$FORESTNAME, cex=0.75, pos=3, col="white")
33
Sample Dataset
Exercise
## Exercise 2
## 2.1 Get projection of nfspoly.
## 2.2 Display the area of interest boundary (bndpoly) with no color fill.
## 2.3 Add the wilderness boundary (wildpoly) with green fill and with no outline.
## 2.4 Get total sum of area of nfspoly.
34
SpatialPixels
35
Sample Dataset
SpatialPixels
# National Elevation Dataset (NED) in meters
elevfn <- "SpatialData/uintaN_elevm.img"
# Elevation (m)
36
Sample Dataset
SpatialPixels
# Forest/Nonforest map
fnffn <- "SpatialData/uintaN_fnf.img"
37
Sample Dataset
SpatialPixels - Importing/Exporting
## Importing SpatialPixels (raster) files
help(raster)
elev <- raster("SpatialData/uintaN_elevm.img")
elev
dim(elev)
class(elev)
## Importing a raster stack
help(stack)
rstack <- stack("SpatialData/uintaN_elevm.img", "SpatialData/uintaN_fnf.img")
rstack
dim(rstack)
class(rstack)
## Exporting a raster
help(writeRaster)
writeRaster(rstack, filename="Outfolder/rstack1.img")
# Add integer datatype and notice size of file
writeRaster(rstack, filename="Outfolder/rstack2.img", datatype="INT1U",
overwrite=TRUE)
38
Sample Dataset
SpatialPixels - Spatial Info
## Get spatial information for elev
elev
class(elev)
inMemory(elev)
unique(elev)
extent(elev)
boxplot(elev)
dim(elev)
res(elev)
ncell(elev)
maxValue(elev)
freq(elev)
#
#
#
#
#
#
#
#
#
Checks if raster is stored in memory (RAM)
Get unique values of raster
Get extent of raster
Display boxplot of raster values
Get dimensions of raster
Get the resolution of raster (x y)
Number of cells of raster
Get Maximum value in raster
Get frequency table of raster
## Overview of functions in raster package
http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/raster/html/raster-package.html
39
Sample Dataset
SpatialPixels - slots
## Get slot names for class SpatialPixels
getClass("SpatialPixels")
slotNames(elev)
str(elev)
# @file
# @data
# @legend
# @title
# @extent
# @ncols
# @nrows
# @crs
the source (file name) of raster
list of attributes of raster
the legend features of raster, if exist
the title of raster, if exists
the bounding box of raster
number of columns of raster
number of rows of raster
coordinate reference system slot
## Get slot names for class SpatialPixels slots
slotNames([email protected])
[email protected]@name
[email protected]@nbands
slotNames([email protected])
[email protected]@min
40
Sample Dataset
SpatialPixels - Displaying
## Display rasters
plot(elev)
# Display with terrain colors – 4 classes
plot(elev, col=terrain.colors(n=4))
# Display with terrain colors – 20 classes
plot(elev, col=terrain.colors(n=20))
# Exclude axis labels
plot(elev, col=terrain.colors(n=4), axes=FALSE)
# Get help with other colors for continuous data
help(terrain.colors)
# Display elevation data with contours
contour(elev)
41
Sample Dataset
SpatialPixels - Deriving
## Derive new rasters
# Generate slope and aspect and hillshade from elevation dataset
slp <- terrain(elev, opt=c('slope'), unit='radians')
asp <- terrain(elev, opt=c('aspect'), unit='radians')
hill <- hillShade(slope, asp)
# Generate hillshade
# Display hillshade in grey scale
plot(hill, col=grey(0:100/100), legend=FALSE)
# Classify raster
elevcl <- cut(elev, breaks=5)
# 5 classes
plot(elevcl)
# Display classified raster
cols <- c("brown", "palegreen", "orange1", "dark green", "snow")
plot(elevcl, col=cols, breaks=c(0:5))
# Converts raster from meters to feet
elevft <- elev * 3.28084
summary(elev)
summary(elevft)
42
Sample Dataset
SpatialPixels - fnf
## Get spatial information for fnf (thematic raster)
fnf <- raster("SpatialData/uintaN_fnf.img")
unique(fnf)
# Get unique values of raster
extent(fnf)
# Get extent of raster
barplot(fnf)
# Display barplot of raster values
## Display fnf
fnfvals <- sort(unique(fnf))
plot(fnf, breaks=fnfvals, col=c("green", "brown", "blue"),
legend=TRUE)
## Overview of functions in raster package
http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/raster/html/raster-package.html
43
Sample Dataset
reclass raster
## Reclass raster layer to 2 categories
help(reclassify)
fromvect <- c(0,1,2,3)
tovect <- c(2,1,2,2)
rclmat <- matrix(c(fromvect, tovect), 4, 2)
fnfrcl <- reclassify(x=fnf, rclmat)
fnfrcl
unique(fnfrcl)
freq(fnfrcl)
## To save to file
help(writeRaster)
stratrcl <- reclassify(x=fnf, rclmat, datatype='INT2U',
filename="SpatialData/uintaN_fnfrcl.img", overwrite=TRUE)
## Plot new reclassed raster
fnfrclvals <- c(0, unique(fnfrcl))
plot(fnfrcl, breaks=fnfrclvals, col=c("green","brown"), legend=TRUE)
barplot(fnfrcl)
44
Temporary Files
## Note:
Functions in raster package create temporary files if the values of an output
RasterLayer cannot be stored in memory (RAM). This can happen when no
filename is provided to a function an in functions where you cannot provide
a filename. Temporary files are automatically removed at the start of each session.
## During session:
showTmpFiles()
removeTmpFiles()
## You can change where the temporary folder is in the Rprofile.site file.
## (C:\Program Files\R\R-3.0.0\etc\Rprofile.site)
# Added
options(scipen=6)
options(rasterTmpDir='c:/Temp/')
rasterOptions()
rasterOptions(chunksize = 1e+04, maxmemory = 1e+06)
45
Sample Dataset
Exercise
## Exercise 3
## 3.1 Get the minimum value of elev
## 3.2 Display elev raster with heat colors and 10 classes
## 3.3 Generate slope (in degrees)
## 3.4 Classify the slope raster from 3.3 into 4 classes
## 3.5 Display the classified slope raster with unique colors. Use colors() to select
color choices. Exclude axes labels.
46
Spatial Help Links
# Summary of tools for Spatial data analysis
http://cran.r-project.org/web/views/Spatial.html
# Presentation – Introduction to representing spatial objects in R – Roger Bivand
http://geostat-course.org/system/files/monday_slides.pdf
# List of common projections
http://www.remotesensing.org/geotiff/proj_list/
# List of common projections
http://www.remotesensing.org/geotiff/proj_list/
## Help for point graphics
http://127.0.0.1:11789/library/graphics/html/points.html
## Overview of functions in raster package
http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/raster/html/raster-package.html
47

similar documents