Land cover (25 m): reclassify categories

Load packages, read data and source custom scripts

rm(list = ls())
library(stars)

#> Loading required package: abind

#> Loading required package: sf

#> Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0

path_proj <- day2day::git_path()
path_data <- file.path(path_proj, "data")
path_raw <- file.path(path_data, "raw")
path_cleaned <- file.path(path_data, "cleaned")

path_lcover <- file.path(path_raw, "land-cover")

lclass_file <- file.path(path_lcover, "landcover_classification_speed_v2020.csv")
lclass <- readr::read_csv(lclass_file)

#> Parsed with column specification:
#> cols(
#>   speed_lc_class = col_character(),
#>   lc_code = col_double(),
#>   source_no = col_double(),
#>   source = col_character(),
#>   land_cover = col_character(),
#>   lc_code_int = col_double(),
#>   speed_lc_name = col_character(),
#>   speed_lc_numeric = col_double(),
#>   ErickClass = col_character(),
#>   ErickCode = col_double()
#> )

lcover_file <- file.path(path_lcover, "landcover_composite_map.tif")
lcover <- read_stars(lcover_file, proxy = TRUE)

Small change on coordinates

Round coordinates to be compatible with other datasets and to make aggregation easier and faster.

dd <- st_dimensions(lcover)
dd$x$offset <- round(dd$x$offset)
dd$y$offset <- round(dd$y$offset)
attr(lcover, "dimensions") <- dd

Define new classes

We create an auxiliary vector where:

• the position of the vector represent the 1000 * lc_code
• the value of the element represent the new code to apply

Using this auxiliary variable make the re-classification faster.

lc_code <- round(1000 * lclass$lc_code)
classes_new <- rep(NA, max(lc_code))
classes_new[lc_code] <- lclass$ErickCode

Reclassify and create new raster

The original land categories from landcover_composite_map.tif does not match exactly with lc_code due to number precision. To match perfectly the integer counterpart in thousands, it is necessary to round the numbers after multiplying by 1000.

st_reclass <- function(X) {
    X[[1]] <- array(classes_new[round(1000 * X[[1]])] , dim(X[[1]]))
    X
}

lcover_new <- stars:::collect(lcover, match.call(), "st_reclass", "X")
system.time(
write_stars(lcover_new, file.path(path_cleaned, "landcover_reclass.tif"), type =
            "Byte")
)

#> ====================================================================================================

#>    user  system elapsed 
#> 177.612  95.113 272.931

Visualize reclassified land cover

lcover_new <- read_stars(file.path(path_cleaned, "landcover_reclass.tif"), proxy = TRUE)

cols <- c("white", "chocolate", "royalblue", "yellow", "forestgreen", "greenyellow", "black", "red")
plot(lcover_new, reset = FALSE, col = cols, main = "Land cover: reclassified")

Time to execute the task

Only useful when executed with Rscript.

proc.time()

#>    user  system elapsed 
#> 196.771  96.555 293.208