Raster overlay on Leaflet map

rasta.leaf.R

# Expanding on the great code from Milos Popovic
#   new display as layer on Leaflet map, instead of ggplot

# original: https://github.com/milos-agathon/mapping-raster-files-with-terra-in-r
################################################################################
#                 Mapping OSM and satelitte data with terra in R
#                 Milos Popovic
#                 2022/09/22
################################################################################

# libraries we need
libs <- c(
    "tidyverse", "sf",
    "osmdata", "terra",
    "httr", "XML", "lwgeom",
    "raster", "leaflet"    # added for leaflet
)

# install missing libraries
installed_libs <- libs %in% rownames(installed.packages())
if (any(installed_libs == F)) {
    install.packages(libs[!installed_libs])
}

# load libraries
invisible(lapply(libs, library, character.only = T))

# Montserrat font
sysfonts::font_add_google("Montserrat", "Montserrat")
showtext::showtext_auto()

# 1. GET DELHI OSM DATA
#----------------------

city <- "Delhi, India"

get_bounding_box <- function() {
    bbox <- osmdata::getbb(city)
    return(bbox)
}

bbox <- get_bounding_box()

road_tags <- c(
    "motorway", "trunk",
    "primary", "secondary",
    "tertiary", "motorway_link",
    "trunk_link", "primary_link",
    "secondary_link", "tertiary_link"
)

get_osm_roads <- function() {
    roads <- bbox %>%
        opq() %>%
        add_osm_feature(
            key = "highway",
            value = road_tags
        ) %>%
        osmdata_sf()

    return(roads)
}

roads <- get_osm_roads()

# plot roads
plot(sf::st_geometry(roads$osm_lines))

# 2. GET BUILT-UP DATA
#----------------------

# website
url <- "https://glad.umd.edu/users/Potapov/GLCLUC2020/Built-up_change_2000_2020/"

get_raster_links <- function() {
    # make http request
    res <- GET(url)
    # parse data to html format
    parse <- XML::htmlParse(res)
    # scrape all the href tags
    links <- XML::xpathSApply(parse, path = "//a", xmlGetAttr, "href")
    # grab links
    lnks <- links[-c(1:5)]
    # make all links and store in a list
    for (l in lnks) {
        rlinks <- paste0(url, lnks)
    }

    return(rlinks)
}

rlinks <- get_raster_links()

# bbox values
delhi_ext <- unname(c(
    bbox[1, ][1], bbox[1, ][2],
    bbox[2, ][1], bbox[2, ][2]
))

delhi_ext

# create Delhi extent
de <- c(77.06194, 77.38194, 28.49172, 28.81172)

get_builtup_data <- function() {
    l <- rlinks[grepl("30N_070E", rlinks)]
    ras <- terra::rast(l)
    delhi_ras <- terra::crop(ras, de)
    df <- terra::as.data.frame(delhi_ras, xy = T)
    names(df)[3] <- "value"
    # define categorical values
    df$cat <- round(df$value, 0)
    df$cat <- factor(df$cat,
        labels = c("no built-up", "new", "existing")
    )

    return(df)
}

df <- get_builtup_data()


# 3. DISPLAY ON LEAFLET MAP (instead of ggplot)
#-----------------------------------------------

# build raster data from matrix
tmp <- t(matrix(df$value, ncol=1280))      # rows are 1280 points long
r <- raster(
  tmp, 
  xmn= range(df$x)[1], xmx= range(df$x)[2],
  ymn= range(df$y)[1], ymx= range(df$y)[2],
  crs="+proj=longlat +datum=WGS84"
)
# plot(r)    # admire the data (optional)

# leaflet does not know 'gray20', changed to '#333333'
colrs <- c("#333333", "#FCDD0F", "#287DFC")   

pal <- leaflet::colorNumeric(colrs, c(0,1,2), na.color= "transparent")
# previewColors(pal, df[1:5,3])

leaflet() |>
  addTiles(group= "OSM (default)") |>
  addProviderTiles("Esri.WorldImagery", group= "Satellite") |>
  addProviderTiles("OpenTopoMap", group= "Topo") |>
  addRasterImage(r, colors = pal, opacity = 0.5, group= 'Construction') |>
  addLegend(colors= colrs, labels= c('not built', 'new', 'existing'),
            title= "Construction<br> 2000-2020") |>
  addLayersControl(
    baseGroups= c("OSM (default)","Satellite","Topo"),
    overlayGroups=c('Construction')
)

Advantages over ggplot include interactivity like zoom, pan, swap layers, opacity control and more.