3D plot using Rayshader

Question

My aim is to create a crisp spike map (e.g. Example) based on the distribution of points in a given country. Although I got the "spikes", I do not know how to include the country's borders. Thanks in advance

My data:

Sample data: (Points)

dput(points_sf)
new("SpatialPoints", coords = structure(c(118.695, 121.97279, 
122.871302, 121.15825, 125.822884, 121.04857, 125.700066, 123.369778, 
120.158337, 120.9861, 120.97923, 120.414345, 121.234863, 124.02077, 
121.04736, 121.10176, 120.100426, 125.04301, 123.169, 121.16667, 
118.77164, 121.30082, 120.52274, 124.10022, 121.07869, 120.21649, 
125.58032, 121.07088, 123.169, 120.76037, 121.068504, 118.75597, 
124.96605, 120.759026, 120.88028, 121.04468, 123.64, 124.144, 
118.66307, 125.574646, 124.547777, 125.38587, 121.05398, 120.95626, 
121.06513, 121.00975, 121.03606, 120.5966, 118.73818, 118.733299, 
9.798928, 11.824005, 9.83088, 14.594145, 7.458123, 14.521641, 
7.293635, 9.485694, 14.973889, 14.5186, 14.58382, 12.659684, 
16.837931, 12.977237, 13.11158, 14.672202, 16.124355, 7.572968, 
8.389, 17.08333, 9.731328, 17.13036, 18.177658, 9.697855, 14.635825, 
15.955174, 7.10541, 14.649635, 8.389, 14.264466, 14.642514, 9.776954, 
8.136726, 14.97685, 15.074081, 14.547942, 13.381, 9.85, 9.745967, 
7.028661, 6.302777, 7.02564, 12.360868, 14.068985, 14.654883, 
14.542219, 14.438277, 16.460638, 9.762018, 9.7333), dim = c(50L, 
2L), dimnames = list(NULL, c("coords.x1", "coords.x2"))), bbox = structure(c(118.66307, 
6.302777, 125.822884, 18.177658), dim = c(2L, 2L), dimnames = list(
 c("coords.x1", "coords.x2"), c("min", "max"))), proj4string = new("CRS", 
 projargs = "+proj=longlat +datum=WGS84 +no_defs"))

Countries shapefile:

countries<- geodata::gadm(
 country = c("PH"),level = 1, path = getwd() ) %>%
 sf::st_as_sf() %>%
 sf::st_cast("MULTIPOLYGON") %>% st_transform(crs = st_crs(points_sf))
countries_geom <- st_geometry(countries)

This is the code I am using:

points_main <- as_Spatial(points_sf$geometry)
point_extent <- extent(points_main) 
buffer <- 0.01 
density_raster <- raster(
 xmn = xmin(point_extent) - buffer, 
 xmx = xmax(point_extent) + buffer, 
 ymn = ymin(point_extent) - buffer, 
 ymx = ymax(point_extent) + buffer, 
 resolution = 0.01, # Set finer resolution to capture point density
 crs = "+proj=longlat +datum=WGS84"
)
density_raster <- rasterize(
 points_main, 
 density_raster, 
 fun = "count", background = 0)
density_raster[density_raster == 0] <- NA
pop_mat <- density_raster |>
 as("Raster") |>
 rayshader::raster_to_matrix()
pop_mat |>
 rayshader::height_shade(texture = texture) |>
 add_overlay(polygon_layer) %>%
 rayshader::plot_3d(
 heightmap = pop_mat,
 solid = F,
 soliddepth = 0,
 zscale = 15,
 shadowdepth = 0,
 shadow_darkness = .95,
 windowsize = c(800, 800),
 phi = 65,
 zoom = .65,
 theta = -30,
 background = "white"
 )