giscoR

install.packages("giscoR")

library(remotes)
install_github("rOpenGov/giscoR")

# Enable this universe
options(repos = c(
  ropengov = "https://ropengov.r-universe.dev",
  CRAN = "https://cloud.r-project.org"
))


install.packages("giscoR")

library(giscoR)
library(sf)
library(dplyr)

# Different resolutions
DNK_res60 <- gisco_get_countries(resolution = "60", country = "DNK") %>%
  mutate(res = "60M")
DNK_res20 <-
  gisco_get_countries(resolution = "20", country = "DNK") %>%
  mutate(res = "20M")
DNK_res10 <-
  gisco_get_countries(resolution = "10", country = "DNK") %>%
  mutate(res = "10M")
DNK_res03 <-
  gisco_get_countries(resolution = "03", country = "DNK") %>%
  mutate(res = "03M")


DNK_all <- bind_rows(DNK_res60, DNK_res20, DNK_res10, DNK_res03)

# Plot ggplot2

library(ggplot2)

ggplot(DNK_all) +
  geom_sf(fill = "tomato") +
  facet_wrap(vars(res)) +
  theme_minimal()

# Labels and Lines available

labs <- gisco_get_countries(
  spatialtype = "LB",
  region = "Africa",
  epsg = "3857"
)

coast <- gisco_get_countries(
  spatialtype = "COASTL",
  epsg = "3857"
)

# For zooming
afr_bbox <- st_bbox(labs)

ggplot(coast) +
  geom_sf(col = "deepskyblue4", size = 3) +
  geom_sf(data = labs, fill = "springgreen4", col = "darkgoldenrod1", size = 5, shape = 21) +
  coord_sf(
    xlim = afr_bbox[c("xmin", "xmax")],
    ylim = afr_bbox[c("ymin", "ymax")]
  )

ITA <- gisco_get_nuts(country = "Italy", nuts_level = 1)

ggplot(ITA) +
  geom_sf() +
  geom_sf_text(aes(label = NAME_LATN)) +
  theme(axis.title = element_blank())

# Get shapes
nuts3 <- gisco_get_nuts(
  year = "2016",
  epsg = "3035",
  resolution = "3",
  nuts_level = "3"
)

# Group by NUTS by country and convert to lines
country_lines <- nuts3 %>%
  group_by(
    CNTR_CODE
  ) %>%
  summarise(n = n()) %>%
  st_cast("MULTILINESTRING")


# Use eurostat
library(eurostat)

popdens <- get_eurostat("demo_r_d3dens")
popdens <- popdens[popdens$time == "2018-01-01", ]


# Merge data
nuts3.sf <- merge(nuts3,
  popdens,
  by.x = "NUTS_ID",
  by.y = "geo",
  all.x = TRUE
)

# Breaks and labels

br <- c(0, 25, 50, 100, 200, 500, 1000, 2500, 5000, 10000, 30000)

nuts3.sf$values_cut <- cut(nuts3.sf$values,
  breaks = br,
  dig.lab = 5
)

labs_plot <- prettyNum(br[-1], big.mark = ",")


# Palette
pal <- hcl.colors(length(br) - 1, "Lajolla")


# Plot

ggplot(nuts3.sf) +
  geom_sf(aes(fill = values_cut), size = 0, color = NA, alpha = 0.9) +
  geom_sf(data = country_lines, col = "black", size = 0.1) +
  # Center in Europe: EPSG 3035
  coord_sf(
    xlim = c(2377294, 7453440),
    ylim = c(1313597, 5628510)
  ) +
  labs(
    title = "Population density in 2018",
    subtitle = "NUTS-3 level",
    caption = paste0(
      "Source: Eurostat, ", gisco_attributions(),
      "\nBased on Milos Popovic: https://milospopovic.net/how-to-make-choropleth-map-in-r/"
    )
  ) +
  scale_fill_manual(
    name = "people per sq. kilometer",
    values = pal,
    labels = labs_plot,
    drop = FALSE,
    guide = guide_legend(
      direction = "horizontal",
      keyheight = 0.5,
      keywidth = 2.5,
      title.position = "top",
      title.hjust = 0.5,
      label.hjust = .5,
      nrow = 1,
      byrow = TRUE,
      reverse = FALSE,
      label.position = "bottom"
    )
  ) +
  theme_void() +
  # Theme
  theme(
    plot.title = element_text(
      size = 20, color = pal[length(pal) - 1],
      hjust = 0.5, vjust = -6
    ),
    plot.subtitle = element_text(
      size = 14,
      color = pal[length(pal) - 1],
      hjust = 0.5, vjust = -10, face = "bold"
    ),
    plot.caption = element_text(
      size = 9, color = "grey60",
      hjust = 0.5, vjust = 0,
      margin = margin(t = 5, b = 10)
    ),
    legend.text = element_text(
      size = 10,
      color = "grey20"
    ),
    legend.title = element_text(
      size = 11,
      color = "grey20"
    ),
    legend.position = "bottom"
  )

gisco_set_cache_dir("./path/to/location")