In arc2r/book: All things R Spatial or A guide to spatial R for Recovering ArcGIS Users

//todo reintegrate this into the document

Adding Raster Maps Files (geotiffs)

Its complicated enough with standard Raster Data, and even more so with raster maps. There are two packages that facilitate the combination of raster maps and ggplot:

• rasterVis
• RStoolbox

However, we will do without these packages to learn more on how rasters work and be as flexible as possible.

If you don't have the time or patience to read through the explanation, the following function takes a path to a raster and returns a dataframe with x/y coordinates as well as a column with the hex colours for the according cell.

rastermap_to_dataframe <- function(path) {
  require(tiff)
  require(raster)
  raster_file <- raster::brick(path)
  bands_n <- nbands(raster_file)
  if(bands_n == 1){
    raster_indexed <- readTIFF(path,indexed = TRUE)
    raster_colormap <- attr(raster_indexed,"color.map")
    raster_colormap_hex <- rgb(raster_colormap[1,],raster_colormap[2,],raster_colormap[3,])
    raster_df <- data.frame(coordinates(raster_file), values = raster_colormap_hex[getValues(raster_file)+1])
  } else if(bands_n == 3){
    raster_df <- data.frame(
      coordinates(raster_file),
      values = rgb(values(raster_file[[1]]),values(raster_file[[2]]),values(raster_file[[3]]),maxColorValue = 255)
    )
  }
  return(raster_df)
}

With Colours

Tif files with colours can either be single band or include multiple bands in a single tif file. These need to be treated differently. Let's take two variants of the swiss raster map 1:50k map. One is a single band dataset, the other a 3-band dataset.

library(raster)
library(rgdal)

path_1b <- "sample_data/SMR_Musterdaten/SMR50_LV95_KOMB_Mosaic.tif"
path_3b <- "sample_data/SMR_Musterdaten/SMR50_LV95_KREL_10L_Mosaic.tif"

raster(path_1b)
raster(path_3b)

As you can see in the output of raster(path_3b), only one of the 3 bands was imported with raster::raster. In the case of multiband rasters, it's better to use raster::brick:

map_1b <- raster(path_1b)
map_3b <- brick(path_3b)

Singleband

A singleband RasterLayer can be regarded as vector of numeric values arranged in a matrix with some additional metadata. The metadata can be view with str(), the values can be obtained with getValues().

str(map_1b)
head(getValues(map_1b),50)

extent(map_1b)

The first 50 cells of this RasterLayer contain the value 37. This value 37 is associated with a specific colour, which can be viewed in the colortable of the RasterLayer.

coltab <- colortable(map_1b)
head(coltab)
scales::show_col(coltab,labels = FALSE)

In order to plot this RasterLayer in ggplot, we need to convert it into a dataframe. You can imagine that this is converting the wide matrix into a very long table with three columns: x and y coordinates of the lower left corner of each cell and a value specifying what the cell holds.

map_1b_df <- data.frame(coordinates(map_1b),values = getValues(map_1b))
head(map_1b_df)

We can now use the values column in ggplot2 to specify the fill aesthetic. If we then pass the comfortable as our fill values, the colours will be matched correctly. However we need two additional tricks:

• name the colourtable to make sure the correct colours are correct to the according values
• wrap the values column in factor() so that the ggplot2 does not mistake the numeric values as a continuous scale:

names(coltab) <- 0:255 # 1:256 would be wrong by 1

ggplot() +
  geom_raster(data = map_1b_df, aes(x,y, fill = factor(values))) +
  scale_fill_manual(values = coltab) +
  theme(legend.position = "none") +
  coord_equal()

Note that there are thee ways to plot this type of data in ggplot2: geom_rect, geom_tile and geom_raster. The latter is "a high performance special case for when all the tiles are the same size" (see the docs.

Alternatively, you can add the hex colourcode to the dataframe and use scale_fill_identitiy().

map_1b_df <- mutate(map_1b_df, colorvalue = coltab[values+1])
ggplot() +
  geom_raster(data = map_1b_df, aes(x,y, fill = colorvalue)) +
  scale_fill_identity() +
  theme(legend.position = "none") +
  coord_equal()

The big drawback of using ggplot2 in this way: you cannot add a second layer with a fill aesthetic, since ggplot2 only allows one per plot.

Multiple bands

In case of a 3-Band Raster Brick, each layer holds a value of 0 to 255 representing colours in the RGB colour model. Usually, the first layer represents Red, the second Green and the third Blue. We can subset the RasterBrick into the individual RasterLayers using double brackets ([[1]] for the first layer).

map_3b

map_3b[[1]]

As with the singleband raster, we can extract the Values (this time RGB-Values) of each layer with getValues().

head(values(map_3b[[1]]),50)

Again like with the single band raster, we now transform the Raster into a dataframe, this time storing the rgb values in columns (instead of the hex values as before).

map_3b_df <- data.frame(
  coordinates(map_3b),
  r = values(map_3b[[1]]), 
  g = values(map_3b[[2]]),
  b = values(map_3b[[3]]))

head(map_3b_df)

With the function rbg() we can transform values into hex-codes.

rgb(10,10,10,maxColorValue = 255)

ggplot() + 
  geom_raster(data = map_3b_df, aes(x,y,fill = rgb(r,g,b,maxColorValue = 255))) + 
  scale_fill_identity() +
  coord_equal()

Tip: if your file is to large and ggplot has trouble plotting it fast, you can use aggregate to reduce the resolution

map_3b_lowres <- raster::aggregate(map_3b, fact = 5)

Greyscale

In greyscale maps, raster cannot find a associated colourtable to match the values to rgb or hex colours.

map_grey_path <- "sample_data/SMR_Musterdaten/SMR50_LV95_KGRS_Mosaic.tif"
map_grey <- raster(map_grey_path)

head(getValues(map_grey))
head(colortable(map_grey))

However, the package tiff can find the rgb values of a tiff, we just need to set indexed = TRUE and get the color.map attribute of the output:

library(tiff)

index <- readTIFF(map_grey_path,indexed = TRUE)
colormap <- attr(index,"color.map")


colormap[,1:8] # only showing the first 8 columns (of 256)

This is an index with three rows (Red, Green and Blue) and 265 columns (with are the Values from getValues()). We can now use rgb() to turn this matrix into a vector of 265 hex colours:

coltab <- rgb(colormap[1,],colormap[2,],colormap[3,])

scales::show_col(coltab,labels = FALSE)

Like before, we now turn the RasterLayer into a dataframe and use the Values and the newly created colourtable to assign the correct values:

map_grey_df <- data.frame(coordinates(map_grey),values = getValues(map_grey))


names(coltab) <- 0:255

ggplot() +
  geom_raster(data = map_grey_df, aes(x,y, fill = factor(values))) +
  scale_fill_manual(values = coltab) +
  theme(legend.position = "none") +
  coord_equal()

arc2r/book documentation built on March 5, 2021, 2:10 p.m.