diff --git a/trajectory_analysis_full.Rmd b/trajectory_analysis_full.Rmd
index 73915bd..8d051b6 100644
--- a/trajectory_analysis_full.Rmd
+++ b/trajectory_analysis_full.Rmd
@@ -940,6 +940,14 @@ library(sf)
 
 ```{r}
 get_shp <- function(ISO2) {
+  # Obtener shp_1
+  shp_0 <- geodata::gadm(
+    ISO2,
+    path = "01_data/raw/shp/",
+    level = 0
+  ) %>%
+    st_as_sf()
+  
   # Obtener shp_1
   shp_1 <- geodata::gadm(
     ISO2,
@@ -969,6 +977,7 @@ get_shp <- function(ISO2) {
 
   # Devolver la lista con shp_1 y shp_2
   list(
+    shp_0 = shp_0,
     shp_1 = shp_1,
     shp_2 = shp_2
   )
@@ -976,6 +985,11 @@ get_shp <- function(ISO2) {
 
 
 shp_list <- unique(gbmt_fit_final$gbmt_fit_total[[1]]$ISO2) %>%
+  append(
+    c("GUY", "SUR", "BOL", "PRY",
+      "URY", "ECU", "BLZ", "HND",
+      "TTO", "VEN")
+  ) %>% 
   set_names() %>%
   map(get_shp)
 ```
@@ -1182,4 +1196,203 @@ walk2(
   x_names = gbmt_fit_final$x_names,
   scale_name = gbmt_fit_final$scale_name
 )
+```
+
+
+# Mapa con puntos
+
+
+```{r}
+join_data_shp2 <- function(data, shp_list) {
+  # Proceso de unión con shapefiles
+  data %>%
+    mutate(
+      pubsalid1 = as.numeric(as.character(pubsalid1)),
+      # City = str_to_lower(City),
+      group = factor(group)
+    ) %>%
+    group_nest(ISO2, country) %>%
+    drop_na(ISO2) %>%
+    mutate(
+      country = str_to_title(country),
+      shp_1 = map(ISO2, ~ shp_list[[.x]]$shp_1),
+      shp_2 = map(ISO2, ~ shp_list[[.x]]$shp_2),
+      shp_inner = map2(
+        data,
+        shp_2,
+        ~ inner_join(
+          .x,
+          .y,
+          by = join_by(pubsalid1)
+        ) %>%
+          st_as_sf()
+      )
+    )
+}
+```
+
+```{r}
+gbmt_interest <- gbmt_fit_final %>%
+  filter(
+    scale_name == "logarithmic",
+    ng == 4
+  ) %>%
+  slice(1)
+
+# Calcula el delta año tras año para cada ciudad
+gbmt_interest$gbmt_fit_total[[1]] <- gbmt_interest$gbmt_fit_total[[1]] %>%
+  arrange(pubsalid1, year) %>%
+  group_by(ISO2, country, group, pubsalid1) %>%
+  mutate(
+    delta_population_imp = population_imp - lag(population_imp),
+    delta_bectuareal1ux_imp = bectuareal1ux_imp - lag(bectuareal1ux_imp),
+    delta_population_imp_norm = population_imp_norm - lag(population_imp_norm),
+    delta_bectuareal1ux_imp_norm = bectuareal1ux_imp_norm - lag(bectuareal1ux_imp_norm)
+  ) %>%
+  summarise(
+    median_delta_population = median(delta_population_imp, na.rm = TRUE),
+    median_delta_bectuareal = median(delta_bectuareal1ux_imp, na.rm = TRUE),
+    median_delta_population_norm = median(delta_population_imp_norm, na.rm = TRUE),
+    median_delta_bectuareal_norm = median(delta_bectuareal1ux_imp_norm, na.rm = TRUE)
+  ) %>%
+  ungroup() %>%
+  mutate(
+    tercile_median_population = ntile(median_delta_population, 3),
+    tercile_median_bectuareal = ntile(median_delta_bectuareal, 3),
+    tercile_median_population_norm = ntile(median_delta_population_norm, 3),
+    tercile_median_bectuareal_norm = ntile(median_delta_bectuareal_norm, 3)
+  )
+
+gbmt_shp_interest <- map(
+  gbmt_interest$gbmt_fit_total,
+  ~ join_data_shp2(.x, shp_list)
+)
+```
+
+```{r}
+country_data <- map(shp_list, 1) %>% 
+  bind_rows()
+
+level1_data <- map(shp_list, 2) %>% 
+  bind_rows()
+```
+
+```{r}
+# country_data <- gbmt_shp_interest[[1]] %>%
+#   select(ISO2, country, shp_1) %>%
+#   unnest(cols = c(shp_1)) %>%
+#   st_as_sf()
+
+city_data <- gbmt_shp_interest[[1]] %>%
+  select(ISO2, country, shp_inner) %>%
+  unnest(cols = shp_inner) %>%
+  st_as_sf() %>%
+  st_centroid()
+```
+
+
+```{r}
+library(biscale)
+```
+
+
+
+```{r}
+city_data_biscale <- bi_class(city_data,
+  x = median_delta_population,
+  y = median_delta_bectuareal,
+  style = "quantile", dim = 3
+)
+
+custom_pal_teal <- c(
+  "1-1" = "#b2d8d8", # teal green for low x, low y
+  "2-1" = "#ba8890",
+  "3-1" = "#9e3547",
+  "1-2" = "#8aa6c2",
+  "2-2" = "#7a6b84",
+  "3-2" = "#682a41",
+  "1-3" = "#4279b0",
+  "2-3" = "#3a4e78",
+  "3-3" = "#311e3b"
+)
+
+
+
+
+ggmap <- level1_data %>%
+  ggplot() +
+  geom_sf(linewidth = 0.05) +
+  geom_sf(
+    data = country_data,
+    linewidth = 0.35, 
+    color = "#3f3f3f",
+    fill = NA
+  ) +
+  geom_sf(
+    data = city_data_biscale,
+    mapping = aes(color = bi_class),
+    #color = "white", 
+    size = 2.5,
+    show.legend = FALSE
+  ) +
+  bi_scale_color(pal = custom_pal_teal, dim = 3) +
+  # bi_scale_color(pal = "GrPink", dim = 3) +
+  bi_theme()
+
+legend <- bi_legend(
+  pal = custom_pal_teal,
+  dim = 3,
+  xlab = "Higher Δ Population",
+  ylab = "Higher Δ Urban Area",
+  size = 6
+)
+
+library(cowplot)
+
+finalPlot <- ggdraw(ggmap) +
+  # draw_plot(ggmap, 0, 0, 1.5, 1.5) +
+  draw_plot(legend, 0.15, 0.25, 0.25, 0.25)
+
+ggsave("ggmap.png",
+       finalPlot,
+       height = 7,
+       width = 7,
+       dpi = 600)
+```
+
+
+```{r}
+ggmap2 <- country_data %>%
+  group_by(country) %>%
+  summarise() %>%
+  ggplot() +
+  geom_sf(fill = "white") +
+  geom_sf(
+    data = city_data_biscale,
+    mapping = aes(color = bi_class),
+    # color = "white",
+    size = 2.5,
+    show.legend = FALSE
+  ) +
+  bi_scale_color(pal = "GrPink", dim = 3) +
+  bi_theme()
+
+legend <- bi_legend(
+  pal = "GrPink",
+  dim = 3,
+  xlab = "Higher Δ Population",
+  ylab = "Higher Δ Urban Area",
+  size = 6
+)
+
+library(cowplot)
+
+finalPlot <- ggdraw() +
+  draw_plot(ggmap, 0, 0, 1, 1) +
+  draw_plot(legend, 0.1, .28, 0.25, 0.25)
+
+ggsave("ggmap.png",
+  finalPlot,
+  dpi = 600
+)
 ```
\ No newline at end of file