Introduction to map making with R
Like basic plots, for making maps there are number packages and functions available. Here we focus on plot and spplot
We will plot raster (elevation) and vector data (administrative boundaries, rivers, lakes and protected areas) of Tanzania. First we will use the plot function. The data are in the zip file on canvas.
Put the answers to the questions in a new RMD file. Include code from this file that you need to write code that runs.
library(sp)
library(raster)
unzip("week4_1_data.zip")
v0 <- readRDS('TZA_adm0.rds')
plot(v0, col= "blue", lwd = 2, border = "red")
title(main = "Tanzania national boundary", sub = "source GADM")
Exercise 1: Fill the polygon using a different color and change increase the thickness of the border
Color individual regions
We will use the region boundaries from GADM database (level = 1)
v2 <- readRDS('TZA_adm2.rds')
library(RColorBrewer)
# Each region should have different color; we will use the color ramp
n <- length(v2$NAME_2)
plot(v2, col=rainbow(n), main = 'Administrative boundary: level 2')
# Now add name of the indiviudal regions to the map
text(v2, v2$NAME_2, cex=0.5)
Exercise 2: Read and plot level 2 administrative boundaries of Tanzania. Give different colors to the regions.
Add lines
Now, let’s add some rivers.
rv <- shapefile('Tanzania_rivers.shp')
plot(v0, border = "blue", lwd = 2, axes = TRUE)
plot(rv, col='pink', add = TRUE)
title(main = "Tanzania rivers")
# Add some more details
legend('bottomleft', legend = 'River', lty = 1, lwd = 2, col = 'blue', bty = "n")
Exercise 3. Change the country boundary color, river color and place the legend in the bottom left of the plot.
Additional polygons
Below we plot the lake and protected areas of Tanzania. First read the required boundaries.
lake <- shapefile('TZA_glwd.shp')
protarea <- shapefile('TZA_wdpa.shp')
plot(v0, lwd = 2, axes = TRUE)
plot(lake, col='lightblue', add = TRUE)
plot(protarea, col='lightgreen', add = TRUE)
title(main = "Tanzania lakes and protected area")
Note the border=“transparent” option suppresses the plotting of polygon borders.
Exercise 4. What happens when you don’t use the border argument in while plotting lake and protected areas?
All the borders will show as default instead of using border=“transparent” which will not show the borders.
Plot Raster
We will plot the elevation of Tanzania.
alt <- raster('TZA_alt.tif')
plot(alt, col = terrain.colors(20), axes=FALSE)
title(main = "Elevation (in m)")
To improve the visualization, we can clamp the higher and lower elevations.
alt <- clamp(alt, 0, 3000)
plot(alt, col = terrain.colors(20), axes=FALSE)
title(main = "Elevation (in m)")
Plotting two rasters
Show elevation and average annual temperature in the same plot.
temp <- raster('TZA_bio1.tif')
temp[temp < 0] = 0
par(mfrow=c(1,2))
plot(alt, col = terrain.colors(20), axes=FALSE)
title(main = "Elevation (in m)")
plot(temp, col = rev(heat.colors(50)), axes=FALSE)
title(main ="Annual Mean Temperature (?C)")
Plotting vectors on top of a raster
We will add the lake, rivers and administrative boundaries to the elevation raster plot
plot(alt, col = terrain.colors(20), legend = FALSE)
plot(lake, col='skyblue1', border = 'transparent', add = TRUE)
plot(rv, col='blue1', add = TRUE)
plot(v0, lwd = 2, axes = TRUE, add = TRUE)
title(main = "Lakes and rivers of Tanzania")
spplot
We will use spplot for different plotting applications.
# Use soil properties information; Soil organic carbon content and soil pH
orc <- raster('TZA_ORC.tif')
ph <- raster('TZA_pH.tif')
soil <- stack(orc, ph)
orc[orc>80]<- 80
spplot(orc, main = list(label="Soil organic carbon content", cExercise = 1))
# change the ph values between 0 to 14
ph <- ph/10
spplot(ph, main = list(label="Soil pH", cExercise = 1))
Now to change the color ramp and change legend position:
brks <- seq(0,60,0.5)
spplot(orc,
at = round(brks, digits=2),
col.regions = rev(terrain.colors(length(brks))), colorkey = list(space = "bottom"),
main = list(label="Soil organic carbon content", cExercise = 1))
brks2 <- seq(4,8,0.5)
spplot(ph,
at = round(brks2, digits=2),
col.regions = rev(terrain.colors(length(brks2))), colorkey = list(space = "bottom"),
main = list(label="Soil pH", cExercise = 1))
Exercise 5. Make similar changes to pH plot. Hint: the brks will take place at different sequence based on pH values
Below we will add the lake and administrative boundary to the elevation raster plot
pols <- list("sp.lines", as(v0, 'SpatialLines'))
brks <- seq(0,60,0.5)
spplot(orc,
sp.layout=pols,
at = round(brks, digits=2),
col.regions = rev(terrain.colors(length(brks))), colorkey = list(space = "bottom"),
main = list(label="Soil organic carbon content", cExercise = 1))
See the difference with plot function. There is no use of add=TRUE argument.
You can also add multiple objects to the plot.
pols1 <- list("sp.lines", as(v0, 'SpatialLines'), col = gray(0.4), lwd = 0.5)
pols2 <- list("sp.polygons", as(lake, 'SpatialPolygons'), fill = 'skyblue1',col="transparent", first = FALSE)
brks <- seq(0,60,0.5)
spplot(orc,
sp.layout=list(pols1, pols2),
at = round(brks, digits=2),
col.regions = rev(terrain.colors(length(brks))), colorkey = list(space = "bottom"),
main = list(label="Soil organic carbon content", cExercise = 1))
pols1 <- list("sp.lines", as(v0, 'SpatialLines'), col = gray(0.4), lwd = 0.5)
pols2 <- list("sp.polygons", as(lake, 'SpatialPolygons'), fill = 'skyblue1',col="transparent", first = FALSE)
pols3 <-list("sp.polygons", as(protarea, 'SpatialPolygons'), col = 'red', first = FALSE)
brks <- seq(0,60,0.5)
spplot(orc,
sp.layout=list(pols1, pols2, pols3),
at = round(brks, digits=2),
col.regions = rev(terrain.colors(length(brks))), colorkey = list(space = "bottom"),
main = list(label="Soil organic carbon content", cExercise = 1))
Exercise 6. Add the protected area boundaries to the plot made above. Do not fill the protected area polygon with any color.
Plot multi-layer raster object
You can also use ssplot to plot multi-layer raster data
brks <- seq(0,70,0.5)
ph <- ph*10
soil <- stack(orc,ph)
spplot(soil,
layout = c(2,1),
at = round(brks, digits=2),
col.regions = rev(terrain.colors(length(brks))), colorkey = list(space = "bottom"),
main = list(label="Soil properties", cExercise = 1))
- Exercise 7.* Save your favorite plot to a “png” file, by filling in the dots.
png(filename = "myplot.png", width = 250, height = 200, units = "mm", res = 300)
spplot(orc,
sp.layout=list(pols1, pols2, pols3),
at = round(brks, digits=2),
col.regions = rev(terrain.colors(length(brks))), colorkey = list(space = "bottom"),
main = list(label="Soil organic carbon content", cExercise = 1))
dev.off()## quartz_off_screen
## 2