Reclassifying Raster Maps in R
Goals
In this exercise we will
- learn to interpret raster map values, in this case representing
landcover types
- reclassify raster map values to simplify our landcover
categories
- add labels to map values to make our final map easier to
interpret
- and plot our map on a basemap using the tmap package
library(readr)
library(dplyr)
library(terra)
library(tmap)
library(tmaptools)
library(maptiles)
library(OpenStreetMap)
library(osmdata)
library(rJava)
library(xlsx)Loading and checking the data
We will be working with crop inventory data from Agriculture and
Agri-food Canada. We have extracted crop inventory data for 10.5 kn
circular buffer around our sampling points near Grande Prairie, Alberta,
Canada.
First let’s import the map and check it’s description
gp2014 <- terra::rast("./data/crop_inventory_2014_GrandePrairie_masked.tif")
gp2014## class : SpatRaster
## dimensions : 9326, 9332, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : -13226340, -12946380, 7253918, 7533698 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
## source : crop_inventory_2014_GrandePrairie_masked.tif
## name : crop_inventory_2014_GrandePrairie_masked
Let’s plot the map. Here we see the circular buffers with a radius of
10.5 km where different values represent different landcover types. Many
of the buffers overlap, which is why we see some blobs instead of
perfect circles, and that’s totally okay.
plot(gp2014)
Let’s now take a look at the attribute table of this map and find out
which values represent which landcover types.
at_gp2014 <- read.csv("./data/attribute_table_crop_inventory_2014_GrandePrairie_masked.csv")
at_gp2014
This map is what I like to call psuedo-categorical. The “Value” column,
which are the values populating the raster, are numerical
representations of the categorical “ClassName” column, which indicate
the landcover type in a meaningful way.
Before we go further, let’s make sure that the map doesn’t contain any
values that don’t have a match in the attribute table.
unique_values_gp2014 <- unique(values(gp2014))
unique_values_at_gp2014 <- unique(at_gp2014$Value)
print(setdiff( unique_values_gp2014, unique_values_at_gp2014 ))## [1] NAprint(setdiff( unique_values_at_gp2014, unique_values_gp2014 ))## integer(0)# Good, no issues!
When we plot the map, it’s a little hard to read as is.
Let’s associate the ClassName variable in our table, to the values in
the raster using the terra::levels function.
To do that, we need to select just the Value and ClassName
variables.
landcover_classes <- at_gp2014 %>%
dplyr::select(Value, ClassName)
Now let’s check to see if there are any levels already in the map.
levels(gp2014)## [[1]]
## [1] ""# Nope. Good, didn't think so!Adding labels
Now, let’s make a copy of our map and add our ClassNames to it.
gp2014_named <- gp2014
levels(gp2014_named) <- landcover_classes
Let’s check the changes . . .
levels(gp2014_named)## [[1]]
## Value ClassName
## 1 2 Spring wheat
## 2 3 Barley
## 3 5 Oats
## 4 14 Canola/rapeseed
## 5 15 Flaxseed
## 6 21 Peas
## 7 41 Grassland
## 8 42 Pasture/forages
## 9 48 Fallow
## 10 49 Too wet to be seeded
## 11 51 Wetland
## 12 52 Water
## 13 53 Shrubland
## 14 55 Coniferous
## 15 56 Broadleaf
## 16 57 Mixedwood
## 17 58 Urban/developed
## 18 59 Exposed land/barren
gp2014_named## class : SpatRaster
## dimensions : 9326, 9332, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : -13226340, -12946380, 7253918, 7533698 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
## source : crop_inventory_2014_GrandePrairie_masked.tif
## categories : ClassName
## name : ClassName
## min value : NA
## max value : NA
Now when we plot the map, we see the ClassNames instead of the values
they represent.
plot(gp2014_named)
Reclassification
Now we can start the reclassification process.
First we need to determine what we want to reclassify.
We have 18 landcover types that we want to simplify into 9 as
follows.
| new_class | new_value | old_class |
|---|---|---|
| Grassland | 100 | Grassland, shurubland, pasture/forages |
| Forest | 101 | Broadleaf, mixedwood, coniferous |
| Barren | 102 | Exposed land/barren, fallow, too wet to be seeded |
| Non-flowering Crop | 103 | Spring wheat, barley, oats, flaxseed |
| Canola | NA | Canola/rapeseed |
| Flowering Crop | NA | Peas |
| Wetland | NA | Wetland |
| Water | NA | Water |
| Developed | NA | Developed/urban |
We only need to make changes to the first 4 categories, the rest can
stay as is.
First we need to make reclassification tables that indicate which values
are getting converted to which. In the table below, values 41, 42, and
53 represent grassland, pasture/forages, and shrubland. We will be
setting all of these values to 100.
grass_rc <- data.frame(from =c(41,42, 53), to =100)
grass_rc
Now we can test the classify function. Let’s save our re-classified map
to a different object name and plot it. If we re-classified the named
version of our map, terra would automatically erase the ClassNames and
create some problems for us, so we use the unnamed map.
gp2014_b <- terra::classify(gp2014, grass_rc)
plot(gp2014_b)
Now that we know how classify works, we can go ahead and do all of our
reclassifications at once, then add the approriate ClassNames to the
reclassified map.
Let’s make the rest of our reclassification tables
forest_rc <- data.frame(from=c(55,56,57), to = 101)
barren_rc <- data.frame(from=c(59,48,49), to = 102)
nonflower_rc <- data.frame(from=c(2,3,5,15), to = 103)
rc_list <- list(grass_rc, forest_rc, barren_rc, nonflower_rc)
rc_data <- dplyr::bind_rows(rc_list)
rc_data
Now let’s reclassify everything
gp2014_rc <- terra::classify(gp2014, rc_data)
gp2014## class : SpatRaster
## dimensions : 9326, 9332, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : -13226340, -12946380, 7253918, 7533698 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
## source : crop_inventory_2014_GrandePrairie_masked.tif
## name : crop_inventory_2014_GrandePrairie_masked
plot(gp2014_rc)
unique_values_gp2014_rc <- unique(values(gp2014_rc))
unique_values_gp2014_rc## crop_inventory_2014_GrandePrairie_masked
## [1,] NaN
## [2,] 14
## [3,] 58
## [4,] 103
## [5,] 101
## [6,] 51
## [7,] 100
## [8,] 102
## [9,] 52
## [10,] 21
Now we can create our new labels
new_landcover_names <- data.frame(values = unique_values_gp2014_rc[-1], Landcover = c("Canola", "Developed", "Non-Flowering_Crop", "Forest", "Wetland", "Grassland", "Barren", "Water", "Flowering_Crop"))
new_landcover_names
and set them to our reclassified map
levels(gp2014_rc)<- new_landcover_names
levels(gp2014_rc)## [[1]]
## values Landcover
## 1 14 Canola
## 2 58 Developed
## 3 103 Non-Flowering_Crop
## 4 101 Forest
## 5 51 Wetland
## 6 100 Grassland
## 7 102 Barren
## 8 52 Water
## 9 21 Flowering_Crop
plot(gp2014_rc)
Let’s quickly make the map in tmap
tmap_mode("view")## ℹ tmap mode set to "view".final_map <-
tm_shape(gp2014_rc) + tm_raster(palette=c("gold3", "grey55","tan4", "olivedrab","lightskyblue3", "darkolivegreen2","darkred", "darkblue","lightpink2"))+
tm_layout(legend.frame = FALSE)+
#tm_compass(position = c("left", "bottom"))+
tm_scale_bar(position = c("left", "bottom"))+
tm_basemap(leaflet::providers$Stadia.Outdoors) # add a basemap so we can see what's around our points more easily##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
## visual variable `col` namely 'palette' (rename to 'values') to col.scale =
## tm_scale(<HERE>).! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.final_map## SpatRaster object downsampled to 3162 by 3164 cells.