Contents

Working with Raster Data


In this tutorial you will learn raster-specific pre-processing and analysis steps such as resampling and reclassification. You will learn how to execute local operations of map algebra in raster calculator and how to summarise raster data using zonal statistics. You will also learn how to create distance, slope and aspect rasters.

We will keep working on siting a solar power station in Zambia, this time using raster data – including the solar power potential.

We will use land cover, and elevation to further optimise the location of the facility: we’d like to avoid cutting trees or building in a swamp. We’d also like to avoid the steepest slopes.

Finally, using solar potential raster (Global Horizontal Irradiance in kWh/m2/year) we will be able to estimate the power which may be available in our selected zones.


Analysis Preparation

Imports

All of these libraries should have been previously installed during the environment set-up, if they have not been installed already you can use install.packages(c("sf", "ggplot2"))

library(sf) # for handling spatial features
library(stars) # for converting between raster and vector data
library(dplyr) # used for data manipulation
library(raster) # useful in some spatial operations
library(ggplot2) # for plotting
library(zeallot) # used for unpacking variables

source('../scripts/helpers.R') # helper script, note that '../' is used to change into the directory above the directory this notebook is in

Linking to GEOS 3.8.1, GDAL 3.1.2, PROJ 7.1.0

Loading required package: abind

Error: package or namespace load failed for ‘stars’ in dyn.load(file, DLLpath = DLLpath, ...):
 unable to load shared object '/srv/conda/envs/notebook/lib/R/library/lwgeom/libs/lwgeom.so':
  libproj.so.15: cannot open shared object file: No such file or directory
Traceback:

1. library(stars)
2. tryCatch({
 .     attr(package, "LibPath") <- which.lib.loc
 .     ns <- loadNamespace(package, lib.loc)
 .     env <- attachNamespace(ns, pos = pos, deps, exclude, include.only)
 . }, error = function(e) {
 .     P <- if (!is.null(cc <- conditionCall(e))) 
 .         paste(" in", deparse(cc)[1L])
 .     else ""
 .     msg <- gettextf("package or namespace load failed for %s%s:\n %s", 
 .         sQuote(package), P, conditionMessage(e))
 .     if (logical.return) 
 .         message(paste("Error:", msg), domain = NA)
 .     else stop(msg, call. = FALSE, domain = NA)
 . })
3. tryCatchList(expr, classes, parentenv, handlers)
4. tryCatchOne(expr, names, parentenv, handlers[[1L]])
5. value[[3L]](cond)
6. stop(msg, call. = FALSE, domain = NA)

Loading Data

We’ll start by again checking to see if we need to download any data

download_data()

Then by reading in a shapefile for Zambia

df_zambia <- read_sf('../data/zambia/zambia.shp')

df_zambia
IDCODECOUNTRYareaQgeometry
761 ZAM Zambia 777795418089 POLYGON ((716452.6 8516873,...

We’ll also load land cover, solar radiance, and elevation raster data from .tif files

solar <- raster('../data/africa/solar.tif')
land_cover <- raster('../data/africa/land_cover.tif')
elevation <- raster('../data/africa/elevation.tif')

plot(solar)
../_images/04) Working with Raster Data_7_0.png

Raster Pre-Processing

In the same way you can clip vectors with a “cookie cutter” outline you can also clip rasters.

df_zambia_4326 <- st_transform(df_zambia, crs=st_crs(solar))
solar_zambia <- crop(solar, df_zambia_4326)

plot(solar_zambia)
plot(df_zambia_4326['geometry'], add=TRUE)
../_images/04) Working with Raster Data_9_0.png 
mask(solar, df_zambia_4326)

Error in sp::CRS(SRS_string = from$wkt): no arguments in initialization list
Traceback:

1. mask(solar, df_zambia_4326)
2. mask(solar, df_zambia_4326)
3. .sf2sp(mask)
4. as(from, "Spatial")
5. asMethod(object)
6. sp::addAttrToGeom(as_Spatial(geom, IDs = row.names(from)), data.frame(from), 
 .     match.ID = FALSE)
7. as_Spatial(geom, IDs = row.names(from))
8. .as_Spatial(from, cast, IDs)
9. sfc2SpatialPolygons(from, IDs)
10. sp::SpatialPolygons(l, proj4string = as(st_crs(from), "CRS"))
11. stopifnot(is(proj4string, "CRS"))
12. is(proj4string, "CRS")
13. as(st_crs(from), "CRS")
14. asMethod(object)
15. CRS_from_crs(from)
16. comment(sp::CRS(SRS_string = from$wkt))
17. sp::CRS(SRS_string = from$wkt)
18. stop(res[[2]])

When we rasterize a Polygon it will try to find an attribute (dataframe column) which can be used to fill in the values of the new raster object, in this case the only numeric column we have relates to the area of a country.

rstr_zambia_4326 <- st_rasterize(df_zambia_4326)

plot(rstr_zambia_4326)
../_images/04) Working with Raster Data_13_0.png
Working with Vector Data Spatial Statistics