Work with GIS data in R

The ducoment is based on the workshop I tought on Oct. 25 2015, Interdepartmental Workshop Series, Oklahoma State University, Stillwater, OK

1. basic

Install and load libraries.

# install.packages("rgdal") 
# install.packages("raster")
# install.packages("dismo")
library("rgdal") # this package is the basis of analyzing GIS data in R; for example, it handle basis coordinate systems, it defines a lot of spatial data types
library("raster")
library("dismo")

2. warm up

R is pretty fun to use, let’s plot some maps from one line of code

require("XML")
# get a map of OK, you could type text that used to search on Google Maps
myMap <- gmap("oklahoma")
raster::plot(myMap)

plot of chunk unnamed-chunk-2

# try different text
plot(gmap("stillwater,OK")) # seems good

plot of chunk unnamed-chunk-3

# get a snapshot of Walmart (satellite image)
plot(gmap("walmart,stillwater,OK",type = "satellite") )

plot of chunk unnamed-chunk-4

# get a snapshot of Boomer Lake (satellite image)
plot(gmap("boomer lake,stillwater,OK",type = "satellite") )

plot of chunk unnamed-chunk-5

3. spatial points

3.1 generate spatial points

# get a map of OK 
myMap <- gmap("Oklahoma",lonlat=TRUE)
# show the extent of this map (raster)
okExtent<- extent(myMap)
okExtent
## class       : Extent 
## xmin        : -105.8468 
## xmax        : -91.56882 
## ymin        : 32.32012 
## ymax        : 38.17342
# generate 10 random points (longitude & latitude) within this extent
myLongitude <- runif(n=10, min=okExtent[1] ,max=okExtent[2] )
myLatitude <- runif(n=10, min=okExtent[3] ,max=okExtent[4] )

# combine longitude and latitude by column
coords <- cbind(myLongitude,myLatitude)
head(coords)
##      myLongitude myLatitude
## [1,]   -93.43081   34.73915
## [2,]   -92.68065   34.33011
## [3,]  -101.66260   34.13352
## [4,]   -95.98800   37.43642
## [5,]   -95.31286   36.34747
## [6,]  -102.96578   33.37099
# make the points spatial
myPoints <- SpatialPoints(coords)

# plot the points
plot(myPoints) # but there is no background

plot of chunk unnamed-chunk-9

# add some background
plot(myMap)
plot(myPoints, add=TRUE, col="red")

plot of chunk unnamed-chunk-10

3.2 generate spatial points with attribute column

# generate random attribute for the points
myAtt <- sample(c("presence","absence"),10,replace=TRUE)

# change myAtt to DataFrame
myAtt <- as.data.frame(myAtt)
head(myAtt)
##      myAtt
## 1  absence
## 2 presence
## 3  absence
## 4  absence
## 5 presence
## 6  absence
# make spatial data frame (spatial points with attributes)
myPoints <- SpatialPointsDataFrame(coords,as.data.frame(myAtt))

# show the attribute of myPoints
myPoints@data
##       myAtt
## 1   absence
## 2  presence
## 3   absence
## 4   absence
## 5  presence
## 6   absence
## 7   absence
## 8   absence
## 9   absence
## 10 presence

3.3 select a subset of points based on attributes

# This is like subsetting a dataframe
myPoints_presence <- myPoints[myPoints$myAtt=="presence", ]

# plot the previously selected points as red
plot(myPoints)
plot(myPoints_presence,add=T,col="red")

plot of chunk unnamed-chunk-12

3.4 save our selected file as a shape file

# the 1st parameter is the object, the 2nd parameter is the path and file name
shapefile(myPoints_presence,filename="temp/output.shp",overwrite=TRUE)
## Error in rgdal::writeOGR(x, filename, layer, driver = "ESRI Shapefile", : Layer creation failed
#file.exists("output.shp")

4. spatial polygons

4.1 build buffer of points, the unit of width depends on the geographic reference system of “myPoints”

# build a dissolved buffer
myBuffer <- buffer(myPoints,width=1)
## Loading required namespace: rgeos
length(myBuffer)
## [1] 1
plot(myBuffer)

plot of chunk unnamed-chunk-14

# build buffer independently for each point
myBuffer <- buffer(myPoints,width=1,dissolve=FALSE)
length(myBuffer)
## [1] 10
plot(myBuffer)

plot of chunk unnamed-chunk-14

# plot them
plot(myMap)

plot of chunk unnamed-chunk-14

plot(myPoints,add=T)
plot(myBuffer,add=T)

plot of chunk unnamed-chunk-14

4.2 load existing polygons (data from http://www.diva-gis.org/Data)

# use the shapefile function, which is used for export if the object is specified
map_state <- shapefile("data/USA_adm1.shp")
## Error: file.exists(extension(x, ".shp")) is not TRUE
# show the summary of the data
summary(map_state)
## Error in summary(map_state): object 'map_state' not found
# show the structure of the data
head(map_state@data, n=5)
## Error in head(map_state@data, n = 5): object 'map_state' not found
# show one colume
map_state$NAME_1
## Error in eval(expr, envir, enclos): object 'map_state' not found
# plot the data
plot(map_state)
## Error in plot(map_state): object 'map_state' not found

4.3 subset

# only select Oklahoma
map_ok <- map_state[map_state$NAME_1 == "Oklahoma", ]
## Error in eval(expr, envir, enclos): object 'map_state' not found
plot(map_ok)
## Error in plot(map_ok): object 'map_ok' not found
# select states large area (> 30)
# step 1: do the logic judgement
selection <- map_state$Shape_Area > 30
## Error in eval(expr, envir, enclos): object 'map_state' not found
# step 2: subset
map_selected <- map_state[selection,]
## Error in eval(expr, envir, enclos): object 'map_state' not found
# the following code shows the same results, but I will nor run it.
#map_selected <- map_state[map_state$Shape_Area < 10,]

# show all columns/fields/attributes of the selections
map_selected@data
## Error in eval(expr, envir, enclos): object 'map_selected' not found
plot(map_selected)
## Error in plot(map_selected): object 'map_selected' not found

4.4 save polygons

shapefile(map_selected,"temp/selected_states.shp",overwrite=TRUE)
## Error in shapefile(map_selected, "temp/selected_states.shp", overwrite = TRUE): object 'map_selected' not found

5. spatial raster

5.1 read/write raster files (data from http://www.worldclim.org)

# read one raster layer
myLayer<- raster("data/bio1.bil")
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
plot(myLayer)
## Error in plot(myLayer): object 'myLayer' not found
# write one raster layer (not run)
# writeRaster(myLayer,filename="temp/ok_bio1.bil",format="EHdr",overwrite=TRUE)

# load several raster layers
# step 1 get a list of file names
list.files("data/") # we need to filter the names
## character(0)
list.files("data/",pattern=".bil") # the names are correct, but we need the full path
## character(0)
list.files("data/",pattern=".bil", full.names = TRUE) # the full name give you a relative path to current working directory
## character(0)
myFiles <- list.files("data/",pattern=".bil", full.names = TRUE) 

# step 2 treat them as raster stack
myLayers <- stack(myFiles)
## Error in x[[1]]: subscript out of bounds
plot(myLayers)
## Error in plot(myLayers): object 'myLayers' not found
# save several raster layers
formattedNames <- paste("temp/",names(myLayers),".bil", sep="")
## Error in paste("temp/", names(myLayers), ".bil", sep = ""): object 'myLayers' not found
formattedNames
## Error in eval(expr, envir, enclos): object 'formattedNames' not found
# not run
# writeRaster(myLayers,filename= formattedNames, format="EHdr", bylayer=TRUE)

5.2 extraction by polygon

# we only want to show Oklahoma, extract raster layer by polygon
raster_ok <- mask(myLayer, map_ok) # we may get error if the reference systems are different
## Error in mask(myLayer, map_ok): object 'myLayer' not found

5.3 projection

# check their CRS
crs(myLayer)
## Error in crs(myLayer): object 'myLayer' not found
crs(map_ok)
## Error in crs(map_ok): object 'map_ok' not found
# unify the CRS
map_ok_new <- spTransform(map_ok, crs(myLayer))
## Error in spTransform(map_ok, crs(myLayer)): object 'map_ok' not found
crs(map_ok_new)
## Error in crs(map_ok_new): object 'map_ok_new' not found
# extract raster by polygon 
raster_ok <- crop( myLayer ,extent(map_ok_new) ) # first cut by a rectangle
## Error in crop(myLayer, extent(map_ok_new)): object 'myLayer' not found
raster_ok <- mask(raster_ok, map_ok_new) # then cut by boundary
## Error in mask(raster_ok, map_ok_new): object 'raster_ok' not found
plot(raster_ok)
## Error in plot(raster_ok): object 'raster_ok' not found

5.4 extract by point

extract(raster_ok,myPoints)
## Error in extract(raster_ok, myPoints): object 'raster_ok' not found

5.5 resample

# we want the new layer to be 3 times coarser at each axis (9 times coarser)
# read current resolution
raster_ok
## Error in eval(expr, envir, enclos): object 'raster_ok' not found
nrow(raster_ok)
## Error in nrow(raster_ok): object 'raster_ok' not found
ncol(raster_ok)
## Error in ncol(raster_ok): object 'raster_ok' not found
extent(raster_ok)
## Error in extent(raster_ok): object 'raster_ok' not found
# define new resolution
newRaster <- raster( nrow= nrow(raster_ok)/3 , ncol= ncol(raster_ok)/3 )
## Error in nrow(raster_ok): object 'raster_ok' not found
# define extent
extent(newRaster) <- extent(raster_ok)
## Error in extent(raster_ok): object 'raster_ok' not found
# fill the new layer with new values
newRaster <- resample(x=raster_ok,y=newRaster,method='bilinear')
## Error in resample(x = raster_ok, y = newRaster, method = "bilinear"): object 'raster_ok' not found
plot(newRaster) # new layer seems coarser
## Error in plot(newRaster): object 'newRaster' not found

5.6 reclassify raster layer

# we want to classify the world into two classes based on temperature, high > mean & low < mean
myLayer<- raster("data/bio1.bil")
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
# values smaller than meanT becomes 1; values larger than meanT will be 2
myMethod <- c(-Inf, 100, 1,  100, Inf, 2)
myLayer_classified <- reclassify(myLayer,rcl= myMethod)
## Error in reclassify(myLayer, rcl = myMethod): object 'myLayer' not found
plot(myLayer_classified)
## Error in plot(myLayer_classified): object 'myLayer_classified' not found

5.7 raster calculation

# read precipitation data 
wet <- raster("data/bio13.bil") # precipitation of wettest month
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
dry <- raster("data/bio14.bil") # precipitation of driest month
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
plot(stack(wet,dry))
## Error in stack(wet, dry): object 'wet' not found
# calculate the difference
diff <- wet - dry
## Error in eval(expr, envir, enclos): object 'wet' not found
#plot(diff)

# calculate the mean of the two month
twoLayers <- stack(wet,dry)
## Error in stack(wet, dry): object 'wet' not found
meanPPT <- calc(twoLayers,fun=mean)
## Error in calc(twoLayers, fun = mean): object 'twoLayers' not found
#plot(meanPPT)

# the following code gives the same results
meanPPT2 <-  (wet + dry)/2
## Error in eval(expr, envir, enclos): object 'wet' not found
# histogram of one layer
hist(twoLayers[[1]])
## Error in hist(twoLayers[[1]]): object 'twoLayers' not found
# correlation between different layers
pairs(twoLayers[[1:2]])
## Error in pairs(twoLayers[[1:2]]): object 'twoLayers' not found

################################

references

  1. Spatial data in R: Using R as a GIS http://pakillo.github.io/R-GIS-tutorial/