1 / 47

Spatial

Spatial. Outline. Types of Spatial Data R Spatial packages S4 R objects Projections Spatial methods SpatialPoints SpatialLines / SpatialPolygons SpatialPixels SpatialPoints - Sample Dataset Importing/Exporting Spatial information and slots Displaying

jana
Download Presentation

Spatial

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Spatial

  2. Outline • Types of Spatial Data • R Spatial packages • S4 R objects • Projections • Spatial methods • 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

  3. 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 spBasic R classes for handling geospatial data maptools Read and write shapefiles. Cannot read the projection file. rgeosInterface to spatial geometry for sp objects rgdalSupports 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

  4. 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)

  5. Spatial Data in R ## Spatial classes (sp package) getClass("Spatial") Bivand, et al. 2008

  6. 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 • bbox – The bounding box, consisting of a matrix of coordinates definingthe extent of the spatial object. • 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

  7. 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"

  8. 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. 2008

  9. SpatialPoints Bivand, et al. 2008

  10. SpatialLines/Polygons List of Objects Object Lists of List Objects Bivand, et al. 2008

  11. SpatialPixels Similar to SpatialPoints with regularly spaced coordinates. Bivand, et al. 2008

  12. Sample Dataset Utah, USA

  13. 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: shrub-montane aspen lodgepole pine spruce-fir alpine

  14. SpatialPoints

  15. 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")

  16. Sample Dataset SpatialPoints – Spatial Info ## Get spatial information for ptshp summary(ptshp) # General info and summary statistics for attribute # data of ptshp mode(ptshp) # Get mode of ptshp class(ptshp) # Get class of ptshp is.projected(ptshp) # Check if ptshp has projection information proj4string(ptshp) # Set or get projection information str(ptshp) # Structure of ptshp

  17. Sample Dataset SpatialPoints - Slots ## Get slot names and access slots for SpatialPoints slotNames(ptshp) # Get name of ptshpcomponents (slots) ptshp@data# Get data frame attributes of ptshp ptshp@coords# Get coordinates of ptshp ptshp@bbox# Get extent of ptshp ptshp@proj4string # Get projection information head(ptshp@data)# Get first 6 rows of data frame of ptshp dim(ptshp) # Get dimensions of data frame of ptshp

  18. 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

  19. 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.

  20. SpatialPolygons

  21. 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"

  22. 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)

  23. 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)

  24. Sample Dataset SpatialPolygons – Spatial Info ## Get spatial information for bndpoly bndpoly # Get general information of polygon summary(bndpoly) # Summary information for polygon mode(bndpoly) # Get mode of bndpoly typeof(bndpoly) # Get typeof class(bndpoly) # Get class of bndpoly proj4string(bndpoly) # Set or get projection information of bndpoly bbox(bndpoly) # Get extent of bndpoly dim(bndpoly) # Get dimension of bndpoly str(bndpoly) # Get structure of bndpoly

  25. Sample Dataset SpatialPolygons - slots (1) (2) (3) # Get slot names for class SpatialPolygons (1) getClass("SpatialPolygons") slotNames(bndpoly) @polygons # the list of Polygons objects @plotOrder# the order to plot the Polygons @bbox# bounding box slot @proj4string # coordinate reference system slot @data # associated data frame (for class SpatialPolygonsDataFrame objects)

  26. 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

  27. Sample Dataset SpatialPolygons – Spatial Info # Get spatial information for nfspoly nfspoly# Get general information of nfspoly summary(nfspoly) # Summary information for nfspoly dim(nfspoly) # Dimensions of nfspoly str(nfspoly) # Get structure of nfspoly

  28. Sample Dataset SpatialPolygons - slots (1) (2) (3) # Get slot names for class SpatialPolygons (1) getClass("SpatialPolygons") slotNames(nfspoly) @polygons # the list of Polygons objects @plotOrder# the order to plot the Polygons @bbox# bounding box slot @proj4string # coordinate reference system slot @data # associated data frame (for class SpatialPolygonsDataFrame objects)

  29. Sample Dataset SpatialPolygons - slots (1) (2) (3) # Get slot names for class Polygons (2) getClass("Polygons") slotNames(nfspoly@polygons[[1]]) # First feature @polygons[[1]]@Polygons # a list of rings (islands/holes) that make up the feature @polygons[[1]]@plotOrder# the order to plot the feature @polygons[[1]]@labpt# label point coordinates of the feature @polygons[[1]]@ID # unique identifier of the feature @polygons[[1]]@area # area of the feature (in units of feature projection)

  30. Sample Dataset SpatialPolygons - slots (1) (2) (3) # Get slot names for class Polygon (3) getClass("Polygon") slotNames(slot(nfspoly@polygons[[1]], "Polygons")[[1]]) @polygons[[1]]@Polygons[[1]] # first ring of feature - class 'Polygon' @polygons[[1]]@Polygons[[1]]@labpt# label point coordinates of the feature ring @polygons[[1]]@Polygons[[1]]@area # area of the feature ring @polygons[[1]]@Polygons[[1]]@coords # coordinates of the feature ring

  31. Sample Dataset SpatialPolygons - slots # Accessing slots of SpatialPolygons nfspoly@data# Data frame attributes of shp nfspoly@polygons# Get info of each polygon within Spatial Polygons nfspoly@plotOrder# Get order of polygons within Spatial Polygons nfspoly@bbox# Get extent of bndpoly # Accessing slots of Polygons slot(nfspoly, "polygons") # (2)slots slotsPolygons <- slot(nfspoly, "polygons") # (2)slots typeof(slotsPolygons) length(slotsPolygons) sapply(slotsPolygons, function(x)slot(x, "ID")) # (2)slots sapply(slotsPolygons, function(x)slot(x, "area")) # (2)slots slotsPolygon <- slot(nfspoly@polygons[[1]], "Polygons") # (3)slots slotsPolygon <- slot(nfspoly@polygons[[2]], "Polygons") # (3)slots

  32. 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)

  33. Sample Dataset SpatialPolygons - Reprojecting ## Notice, the point layer did not display. Check the projections. projection(bndpoly) projection(ptshp) ## Reprojectptshp 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")

  34. 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.

  35. SpatialPixels

  36. Sample Dataset SpatialPixels # National Elevation Dataset (NED) in meters elevfn <- "SpatialData/uintaN_elevm.img" # Elevation (m)

  37. Sample Dataset SpatialPixels # Forest/Nonforest map fnffn <- "SpatialData/uintaN_fnf.img"

  38. 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)

  39. Sample Dataset SpatialPixels - Spatial Info ## Get spatial information for elev elev class(elev) inMemory(elev) # Checks if raster is stored in memory (RAM) unique(elev)# Get unique values of raster extent(elev)# Get extent of raster boxplot(elev) # Display boxplot of raster values dim(elev) # Get dimensions of raster res(elev) # Get the resolution of raster (x y) ncell(elev) # Number of cells of raster maxValue(elev) # Get Maximum value in raster freq(elev) # 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

  40. Sample Dataset SpatialPixels - slots ## Get slot names for class SpatialPixels getClass("SpatialPixels") slotNames(elev) str(elev) # @file the source (file name) of raster # @datalist of attributes of raster # @legend the legend features of raster, if exist # @title the title of raster, if exists # @extent the bounding box of raster # @ncols number of columns of raster # @nrowsnumber of rows of raster # @crs coordinate reference system slot ## Get slot names for class SpatialPixelsslots slotNames(elev@file) elev@file@name elev@file@nbands slotNames(elev@data) elev@data@min

  41. 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)

  42. 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)

  43. 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

  44. 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)

  45. 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)

  46. 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.

  47. 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

More Related