# Map algebra

## Concepts 

The term map algebra typically refers to operations conducted on multiple raster layers using mathematical expressions and boolean logic.  It can be used to build up the complex algorithms used by hydrologists and geomorphologists. Map algebra is also frequently used by ecologists and geographers to build up bespoke calculations of indices that combine many layers, such as human impact, habitat suitability etc. 

## Map algebra in R

The development of the r *raster* package by Robert Hijmans revolutionised the way raster data is handled in R. The raster package makes map algebra simple, intuitive and elegant.  So much so that there is little point in showing any alternative methods in an introductory course, as they are all much more difficult to learn. The only caveat with regards to using the raster package is that although it will handle most jobs it is not well suited to the very largest raster layers. In these cases GRASS remains one of the best open source options for handling operations on extremely large rasters. Postgis can also be used for huge multitemporal raster processing, but this is quite complex and won't be shown on this course.


```{block2, type  ="rmdwarn"}
Always check the size of rasters **before** starting down the path of analysing them. If a raster is too large to fit into computer memory operations in R may hang or crash. The soze of a raster is a combination of the extent and resolution of the layer. The numerical precision of the data also affects the size (integers occupy less memory than floating point numbers). It is  tempting to use the highest resolution data available. This could produce more precise results and more attractive visual output. However do resist the temptation. Use low resolutions at a broad scale and restrict the use of high resolution data to your focal areas. If you really do need to work on very high resolution data at a broad extent you will need to consult someone with experience in lgeoinformatics.
```


```{r, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,message = FALSE,warning = FALSE)
```

### Getting Worldclim layers {#worldclim}



```{r}

require(mapview)
require(leaflet.extras)
library(tmap)
require(raster)
require(rgdal)

mapviewOptions(basemaps =c("OpenStreetMap","Esri.WorldStreetMap", "Esri.WorldImagery", "CartoDB.Positron", "OpenTopoMap", "OpenStreetMap.BlackAndWhite", "OpenStreetMap.HOT","Stamen.TonerLite", "CartoDB.DarkMatter", "Esri.WorldShadedRelief", "Stamen.Terrain"))


```


In addition to providing a simple to use platform for processing raster layers, the raster package (and the dismo package) provide easy access to data sets held on the cloud. The worldclim data set is very commonly used for analysing global climate patterns and is particulary widely used  in species distribution modelling. 

The raster package has an extremely handy and well designed "get_data" function. Not only does this function pull down raster tiles using very inutuitve options to alter the extent and resolution, the function sensibly stores the results locally so that if the file has already been downloaded it is not fetched from the internet again.

The lines below build raster bricks for precipitation, minimum temperature, maximum temperature, elevation (altitude)

```{r}
## Set a low resolution for demo purposes.
library(raster)
path<- "~/geoserver/data_dir/rasters/worldclim"
res<-10
worldclim_prec<-getData(name = "worldclim", var = "prec", res = res,path=path)
worldclim_tmin<-getData(name = "worldclim", var = "tmin", res = res,path = path)
worldclim_tmax<-getData(name = "worldclim", var = "tmax", res = res,path=path)
worldclim_bio<-getData(name = "worldclim", var = "bio", res = res,path=path)
worldclim_alt<-getData(name = "worldclim", var = "alt", res = res,path=path)


```

#### The world altitude layer

We can reuse some code from the terrain analysis section to quickly visualise world elevation at this scale. This only really makes sense for the highest mountains, so the map is centred around the Himalayas.

```{r}
wslope<-terrain(worldclim_alt, opt='slope', unit='radians')
waspect<-terrain(worldclim_alt, opt='aspect',unit='radians')
whill<-hillShade(wslope,waspect,angle=5)

plot(whill, col = gray(0:100 / 100), legend = FALSE, xlim=c(40, 120), ylim=c(10, 60))
plot(worldclim_alt, col = terrain.colors(25), alpha = 0.4, legend = TRUE, add = TRUE, xlim=c(40, 120), ylim=c(10, 60))
```

#### SRTM data for the UK

The raster package also provides an easy way to download and use the shuttle radar telemetry (SRTM) global elevavation layer at 90m resolution. As this coverage has a high resolution at the global extent, the layer is downloaded as tiles. The functon takes a coordinate falling within the tile as input. There are five tiles for the whole of the UK and the Island of Ireland.

```{r,eval=FALSE}
path<- "~/geoserver/data_dir/rasters/srtm"
srtm1<-getData('SRTM', lon=-2, lat=55,path=path)
srtm2<-getData('SRTM', lon=1, lat=55,path = path)
srtm3<-getData('SRTM', lon=-6, lat=55, path= path)
srtm4<-getData('SRTM', lon=-2, lat=60, path = path)
srtm5<-getData('SRTM', lon=-6, lat=60, path = path)

```

The raster package has a useful merge function that joins up the tiles into a continuous tif file. So we can make one like this.

```{r,eval=FALSE}
# srtm<-merge(srtm1,srtm2,srtm3,srtm4,srtm5)
# writeRaster(srtm,file="rasters/srtm_uk.tif")
```

This is a bit slow and the same effective result is also achieved by using a vrt (see next section) 
 
```{r}
path<- "~/geoserver/data_dir/rasters/srtm"
system(sprintf("gdalbuildvrt rasters/srtm.vrt %s/*.tif",path))
```


Now loading the raster layer from a file is very simple.  Just use the raster function to assign it to an object in R. Notice that this is a "lazy" evaluation that does **not** load the raster into memory. The raster declaration simply makes a reference to it. 

Always check on the size of the raster before going any further.

```{r}
# Either
srtm<-raster("rasters/srtm_uk.tif")
# Or 
srtm<-raster("rasters/srtm.vrt")
dim(srtm)

```

The product of the dimensions divided by 1000000 gives the size in MB

```{r}
prod(dim(srtm))/1000000
```

This (i.e `r round(prod(dim(srtm))/1000000,0)` MB) is a manageable file. size in R.  However the actual data is still on disk. There will only be issues with memory if we carry out operations such as trying to place it all into a single mapview.

#### Cropping a raster layer in R {#crop}

Cropping raster layers results in a smaller extent and so lighter demands on memory. The easy way to do this is to use a polygon to set the bounds. Notice that this does not remove pixels falling outside the polygon itself. It crops to the bounding box.

We'll get a polygon for Dorset from the natural earth layer. We don't need to connect to the data base for any of this part.

```{r}
library(rnaturalearth)
library(dplyr)

counties<-ne_states(country = 'United Kingdom', returnclass = c("sf"))
counties %>% filter (name=='Dorset') %>% as("Spatial")->dorset

```


```{r}


dorset_srtm<- raster::crop(srtm,dorset)

dslope<-terrain(dorset_srtm, opt='slope', unit='radians')
daspect<-terrain(dorset_srtm, opt='aspect',unit='radians')
dhill<-hillShade(dslope,daspect,angle=5)

plot(dhill, col = gray(0:100 / 100), legend = FALSE)
plot(dorset_srtm, col = terrain.colors(25), alpha = 0.4, legend = TRUE, add = TRUE) 
plot(dorset,add=TRUE)

```

Let's have a look at this with a mapview.


```{r}
library(RColorBrewer)
pal1<-terrain.colors(100)
pal2<-gray(0:100 / 100)
pal1 <- colorRampPalette(pal1)
pal2 <- colorRampPalette(pal2)

map<-mapview(dorset_srtm,col.regions=pal1, legend=FALSE) + mapview(dhill,col.regions=pal2, legend=FALSE) 
map
```




If we did really want to just work on the SRTM for Dorset we use the mask function to cover areas that are not in Dorset. Any analysis after a mask only uses those pixels.

```{r}
masked_dorset_srtm<-mask(dorset_srtm,dorset)
plot(masked_dorset_srtm)
```

### Raster statistics

The nice part of doing GIS in R is that all the familiar statistical functions are at hand, so any question that you want to ask of the data is just a line of code away. Raster layers are treated just like any other matrix, so we can ask for a summary.


```{r}
summary(masked_dorset_srtm)
```

So the median elevation of the county is 78 meters and the highest point is at 276 meters. Let's find the highlands of Dorset which we'll define as land over 100 meters. This is a simple one liner in R.


```{r}
contours<-rasterToContour(masked_dorset_srtm,levels=c(100))
plot(masked_dorset_srtm>100)
plot(contours, add= TRUE)
```

If we want to know what proportion of Dorset is over 100m we need to know the number of pixels in dorset and the number of pixels over 100m

```{r}
vals<-na.omit(getValues(masked_dorset_srtm)) ## Just take the non masked values

sum(vals>100)/sum(vals>0)
```
So 34% of the land surface of Dorset lies over 100m above sea leval.

Let's try something a little bit harder. Find all the steep(ish) south facing slopes. 

We'll use these definitions.

Steep = slope greater than 5 degrees
South facing = aspect between 130 degrees and 230 degrees

This is a boolean operation. 

```{r}
## Get terrain
slope<-terrain(masked_dorset_srtm, opt='slope', unit='degrees')
aspect<-terrain(masked_dorset_srtm, opt='aspect', unit='degrees')

## Find slopes and set rest to NA
steep_south_slope<-(slope>5 & aspect> 130 & aspect< 230)
steep_south_slope[steep_south_slope<1]<-NA


map + mapview(steep_south_slope,legend=FALSE, na.color ='transparent')
```

Providing layers are of the same extent and the same resolution any sort of mathematical operation can be performed on them together. It is as simple as writing a formula in R to combine a set of vectors.

### Multiple rasters held in a stack (or brick)

A great feature of Hijmans raster package is the use of raster stacks and/or bricks. The idea is similar to a multiband raster, but the elements of a stack or brick can be in separate files. They all must have the same resolution and extent.

```{r}
class (worldclim_prec)

```

The worldclim precipitation object is a raster stack. Raster bricks work in the same way as stacks, but processing a brick can sometimes be wuicker as they are forced into memory. It is usually safer to work with stacks.  You can find the dimenstions with dim.

```{r}
dim(worldclim_prec)
```

Notice that there are three dimenstions. The third is the number of layers, which in the worldclim case is 12, i.e. one for each month of the year. A default print of the stack produces 12 maps.

```{r}
plot(worldclim_prec)
```

The neat feature is that functions work on the stack by driling through the 12 layers. To understand this, let's get a map of total annual rainfall. This is really easy,

```{r}
anprec<-sum(worldclim_prec)
```


```{r}
require(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "BrBG"))
mapview(anprec, col.regions = pal(100), at = c(0, 100, 300,600,1000,3000,10000), legend = TRUE)
```

### Calculating the mean monthly temperature

Two stacks with the same number of layers can also be combined in a really simple way. So a rough mean monthly temperature can be found by adding the maximum and minimum and dividing by 2. The original Worldclim layers hold temperatures as integer values (i.e. the actual value rounded ot one decimal point and multiplied by 10) in orer to save space. When the mean is taken a floating point will result, so at this point we may as well divide by 20 rather than 2 to get the result in degrees centigrade.

```{r}
worldclim_tmean<-(worldclim_tmax + worldclim_tmin) / 20
plot(worldclim_tmean)
```


This is all very neat. It almost looks as if the raster package was designed to make working with worldclim data easy. I wonder why that might be? 


So, mean annual temperature can now be found

```{r}

an_mean_temp<-mean(worldclim_tmean)

```


```{r}
require(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "RdBu")[9:1])
mapview(an_mean_temp, col.regions = pal(100), at = c(-40, 0, 5,10,15,20,25,30), legend = TRUE)
```

## Exercise

The bioclim layers have all been derived from the Worlclim tmax, tmin and prec layers with map algebra. How many of them can you derive?


    BIO1 = Annual Mean Temperature   
    BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp)) 
    BIO3 = Isothermality (BIO2/BIO7) (* 100)   
    BIO4 = Temperature Seasonality (standard deviation *100)  
    BIO5 = Max Temperature of Warmest Month   
    BIO6 = Min Temperature of Coldest Month   
    BIO7 = Temperature Annual Range (BIO5-BIO6)   
    BIO8 = Mean Temperature of Wettest Quarter   
    BIO9 = Mean Temperature of Driest Quarter   
    BIO10 = Mean Temperature of Warmest Quarter    
    BIO11 = Mean Temperature of Coldest Quarter   
    BIO12 = Annual Precipitation   
    BIO13 = Precipitation of Wettest Month   
    BIO14 = Precipitation of Driest Month   
    BIO15 = Precipitation Seasonality (Coefficient of Variation)   
    BIO16 = Precipitation of Wettest Quarter   
    BIO17 = Precipitation of Driest Quarter    
    BIO18 = Precipitation of Warmest Quarter    
    BIO19 = Precipitation of Coldest Quarter  


## Extra bonus: Making climate diagrams for anywhere in the World

One of the motivations behind the worldclim data set is its use in species distribution modelling. The raster extract function adds climate data to any coordinates from WorldClim. So we can add a little utility function to the giscourse package very easily that combines geocoding with Worldclim as a quick look up.

```{r,eval=FALSE}

getwclim<-function(place="Bournemouth",res=2.5) {
path<- "~/geoserver/data_dir/rasters/worldclim"
prec<-getData(name = "worldclim", var = "prec", res = res,path=path)
tmin<-getData(name = "worldclim", var = "tmin", res = res,path = path)
tmax<-getData(name = "worldclim", var = "tmax", res = res,path=path)

require(tmaptools)
pnt<-as(geocode_OSM(place,as.sf=TRUE),"Spatial")
tmins<-unlist(raster::extract(tmin,pnt)/10)
tmaxs<-unlist(raster::extract(tmax,pnt)/10)
precs<-raster::extract(prec,pnt)
df<-data.frame(prec= as.vector(precs),tmax = as.vector(tmaxs),tmin=as.vector(tmins))
df
}

getwclim()
```

```{r,eval=FALSE}
wldiag<-function(df){
 require(climatol)
mat<-rbind(df$prec,df$tmax,df$tmin,df$tmin)
diagwl(mat) 
}


```


```{r}
require(giscourse)
place<-"Bala, Wales"
qmap(place)
getwclim(place) %>% giscourse::wldiag()
```


```{r}
place<-"Mexico city"
qmap(place)
getwclim(place) %>% giscourse::wldiag()
```


```{r,message=FALSE,results='hide',warning=FALSE,echo = FALSE}

rm(list=ls())
lapply(paste('package:',names(sessionInfo()$otherPkgs),sep=""),detach,character.only=TRUE,unload=TRUE)
```


