3. Analysis for Supper.Rmd

---
title: "Story 1: Supper"
output:
  pdf_document: default
  html_document: default
date: "2023-10-09"
---
## 1. Setting up
```{r setup, include=FALSE}
library(GISTools)
library(raster)
library(tmap)
library(rgdal) 
library(maptools) 
library(sp)
library(sf)
library(tidyverse)
library(rgeos)
library(ggplot2)
library(psych)
library(spatstat)
library(spdep)
library(jsonlite)
library(dplyr)
library(gstat)

```

```{r basemap}
subzones <- st_read(dsn = "dataset/basemap", layer = "MP14_SUBZONE_NO_SEA_PL")
#Central area of SG
base_map <- st_read(dsn = "dataset/basemap", layer = "MP14_SUBZONE_NO_SEA_PL")
central_area <- base_map %>% filter(PLN_AREA_N %in% c('DOWNTOWN CORE', 'BUKIT MERAH', 'SINGAPORE RIVER', 'MUSEUM', 'RIVER VALLEY', 'ORCHARD','NEWTON','ROCHOR','TANGLIN', 'KALLANG','OUTRAM')) %>% filter(!SUBZONE_N %in% c( 'CITY TERMINALS','MARITIME SQUARE','ALEXANDRA HILL','TELOK BLANGAH DRIVE','DEPOT ROAD','TELOK BLANGAH WAY','TELOK BLANGAH RISE'))
base_map <- st_transform(central_area, crs=st_crs(central_area)) 
```

### 2. Distribution of Nightlife locations and latenight Restaurants

```{r Nightlife-Restaurants Central distribution}
# Load filtered Nightlife - Restaurants Data

# Load nightlife locations
nightlife_sf <- st_read("dataset/processed_datasets/nightlife_locations.geojson", crs=st_crs(subzones))
nightlife_sf <- st_transform(nightlife_sf, crs=st_crs(subzones)) # project to basemap's CRS

# Load restaurants
restaurants_sf <- st_read("dataset/processed_datasets/restaurants.geojson", crs=st_crs(subzones))
restaurants_sf <- st_transform(restaurants_sf, crs=st_crs(subzones)) # project to basemap's CRS

nightlife_filtered_sf = nightlife_sf
nightlife_filtered_sf$label = 'Nightlife locations'
nightlife_filtered_sf <- nightlife_filtered_sf %>% select(label)

midnight_restaurant_filtered_sf = restaurants_sf %>% filter(midnight_restaurant=="True")
midnight_restaurant_filtered_sf$label = 'Late night restaurants'
midnight_restaurant_filtered_sf <- midnight_restaurant_filtered_sf %>% select(label)

combined <- rbind(nightlife_filtered_sf, midnight_restaurant_filtered_sf)

tm_shape(central_area) + tm_polygons(col="white") + 
  tm_shape(combined) + tm_dots(size= 0.15, col = "label", palette = c("green2", "blue")) + 
  tm_compass(type="8star", size = 2) +
  tm_scale_bar(width = 0.15) +
  tm_layout(legend.format = list(digits = 0),
            legend.position = c("left", "top"),
            legend.text.size = 0.5, 
            legend.title.size = 1,
            title="Nightlife and Latenight restaurants at Central Area",
            title.position = c('left', 'top'))
```

### 3. Hypothesis Testing of 2nd order effect

``` {r PPA}
  
# Point pattern analysis
nightlife_ppp <- as.ppp(nightlife_filtered_sf)
marks(nightlife_ppp) <- "nightlife"
nightlife_ppp$n

latenight_restaurants_ppp <- as.ppp(midnight_restaurant_filtered_sf)
marks(latenight_restaurants_ppp) <- "latenight restaurants"
latenight_restaurants_ppp$n


P <- superimpose(nightlife_ppp, latenight_restaurants_ppp)

# hypothesis testing for the distribution of nightlife and latenight restaurants
nn.sim = vector()
P.r = P
for(i in 1:999){
  marks(P.r) = sample(P$marks)  # Reassign labels at random, point locations don't change
  nn.sim[i] = mean(nncross(P.r[P.r$marks == "nightlife", ] ,P.r[P.r$marks == "latenight restaurants", ])$dist)
}
hist(nn.sim,breaks=30, xlim=c(50,200), main="Histogram of nearest neighbour distance")
actual_dist = mean(nncross(nightlife_ppp,latenight_restaurants_ppp)$dist)
abline(v=actual_dist,col="red")
axis(1, at=actual_dist,labels=round(actual_dist,1))

```

As the p-value is 1, there is strong evidence that the observed distribution of nightlife locations and late-night restaurants is unlikely to have occurred by random chance.

### 4.Crowd Density Interpolation


## Setting up the bar, club and restaurant data frame

```{r bcr data}

base_map <- st_transform(central_area, crs=4326)

#Bars and Clubs dataframe
bcr <- read_csv('dataset/besttime/processed_bars_clubs.csv')
bcr <- na.omit(bcr)
#Restaurants dataframe
res <- read_csv('dataset/besttime/processed_restaurants.csv')
res <- na.omit(res)
res <- res %>% rename(`12am_cd` = `12pm_cd`)

#combining them together
bcr <- bcr %>% select(lat,lon, `10pm_cd`,`11pm_cd`,`12am_cd`,`1am_cd`,`2am_cd`) %>% mutate(cat = "Bars and Clubs")
res <- res %>% select(lat,lon, `10pm_cd`,`11pm_cd`,`12am_cd`,`1am_cd`,`2am_cd`)%>% mutate(cat = "Restaurants")
bcr <- rbind(bcr,res)

#Renaming so that R can recognise the variables
bcr$`cd_10` <- bcr$`10pm_cd`
bcr$`cd_11` <- bcr$`11pm_cd`
bcr$`cd_12` <- bcr$`12am_cd`
bcr$`cd_1` <- bcr$`1am_cd`
bcr$`cd_2` <- bcr$`2am_cd`


#Tranform to spatial dataframe and see how it looks like
coordinates(bcr) <- ~lon + lat

plot(bcr) #notice that 2 points are so much further away than the rest
bcr <- bcr[-c(10,49),] #Remove that one point for better interpolation
plot(bcr) #Corrected plot without annoying points

#used for linear trending
bcr$X <- coordinates(bcr)[,1]
bcr$Y <- coordinates(bcr)[,2]

bcr <- bcr[-zerodist(bcr)[,1],]

#create a spatial layer for basemap
base_spatial <- as(base_map, 'Spatial')
bcr_spatial <- st_as_sf(bcr, coords = c("lon", "lat"))

#quick look at distribution of bars, clubs and restaurants
tmap_mode('plot')
tmap_options(check.and.fix = TRUE)
avail_places <- tm_shape(base_map) + tm_borders(col='black') + tm_fill(col = "white") +
  tm_shape(bcr_spatial) + tm_dots(size = 0.3, col = 'cd_10', title = 'Crowd Density') + # Nightlife and Restaurants locations
  tm_scale_bar(width = 0.15) +
  tm_layout(legend.format = list(digits = 0),
            legend.position = c("left", "top"),
            legend.text.size = 0.75, 
            legend.title.size = 1,
            title="Crowd Density at 10pm",
            title.position = c('left', 'top'))
avail_places  

tmap_save(avail_places, filename = "plots/crowd_density/avail places.png")

```


## Creating raster grid

```{r grid}

#Making the raster grid

grd              <- as.data.frame(spsample(base_spatial, "regular", n=50000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE
fullgrid(grd)    <- TRUE

crs(grd) <- crs(bcr)

```


#### Finding linear trends

So we want to check if any linear trends exist first before we move on to interpolation, to determine whether detrending is needed.

TO be fair, even with little trends present, we will still detrend the data, but this code chunk is also taken as preparation for the next portion.

Lastly, this code chunk is rather long, but we are not surprised considering that we will be plotting for ALL 5 times (10pm,11pm,12am,1am,2am). We wanted to code it into a function, but R doesn't recognise df$x as a vector object :/


```{r trends}

#10pm

#create the linear trend formula
f.bcr <- as.formula(cd_10 ~ X + Y) 

#create the linear trend model based off x and y coordinates
lm.1 <- lm(f.bcr, data = bcr)

#predict the data based off the spatial grid
dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 
r   <- raster(dat.1st)
r.m <- mask(r, base_spatial)

#plot the trend!
l_1 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Linear Trend of Crowd \nDensity at 10pm", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
      tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat',palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants')) +
  tm_legend(legend.outside=TRUE)

#repeat for the rest of the hours
#11pm

f.bcr <- as.formula(cd_11 ~ X + Y) 

lm.1 <- lm(f.bcr, data = bcr)

dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 
r   <- raster(dat.1st)
r.m <- mask(r, base_spatial)

l_2 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Linear Trend of Crowd \nDensity at 11pm", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
        tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants')) +
  tm_legend(legend.outside=TRUE)

#12am

f.bcr <- as.formula(cd_12 ~ X + Y) 

lm.1 <- lm(f.bcr, data = bcr)

dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 
r   <- raster(dat.1st)
r.m <- mask(r, base_spatial)

l_3 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Linear Trend of Crowd \nDensity at 12am", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
        tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants')) +
  tm_legend(legend.outside=TRUE)

#1am

f.bcr <- as.formula(cd_1 ~ X + Y) 

lm.1 <- lm(f.bcr, data = bcr)

dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 
r   <- raster(dat.1st)
r.m <- mask(r, base_spatial)

l_4 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Linear Trend of Crowd \nDensity at 1am", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
        tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants')) +
  tm_legend(legend.outside=TRUE)

#2am

f.bcr <- as.formula(cd_2 ~ X + Y) 

lm.1 <- lm(f.bcr, data = bcr)

dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata=grd))) 
r   <- raster(dat.1st)
r.m <- mask(r, base_spatial)

l_5 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Linear Trend of Crowd \nDensity at 2am", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
        tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants')) +
  tm_legend(legend.outside=TRUE)

l_1
l_2
l_3
l_4
l_5

```

### Creating animations for linear trend


```{r trend animations}


tmap_save(l_1, filename = "plots/crowd_density/linear trend at 10pm.png")
tmap_save(l_2, filename = "plots/crowd_density/linear trend at 11pm.png")
tmap_save(l_3, filename = "plots/crowd_density/linear trend at 12am.png")
tmap_save(l_4, filename = "plots/crowd_density/linear trend at 1am.png")
tmap_save(l_5, filename = "plots/crowd_density/linear trend at 2am.png")

img_list <- image_read(c("plots/crowd_density/linear trend at 10pm.png","plots/crowd_density/linear trend at 11pm.png","plots/crowd_density/linear trend at 12am.png","plots/crowd_density/linear trend at 1am.png","plots/crowd_density/linear trend at 2am.png"))
animation <- image_animate(img_list, fps = 1)  # Adjust fps as needed
image_write(animation, path = "plots/crowd_density/linear trends.gif")


```



## Interpolation for bars and clubs

We repeat the same process for kriging. Similarly, the chunk is very long as we repeat the kriging 5 times :/ We observe that the best model for all would be the wave model, as it seems that the variogram seems to dip after some distance. We manually fit the variogram parameters to each timing to get the best ideal of crowd density at each time.


```{r interpolates}

#10pm
#setting up the formula
f.bcr <- as.formula(cd_10 ~ lon + lat) 

#creating the variogram and fitting it, used wave model
var.smpl.bcr <- variogram(f.bcr, bcr, cloud = FALSE, cutoff=1, width=0.0025)

dat.fit.bcr  <- fit.variogram(var.smpl.bcr, fit.ranges = F, fit.sills = F,
                          vgm(psill=400, model="Wav", range=0.025, nugget=600))

# The following plot allows us to assess the fit
vari_plot_10 <- plot(var.smpl.bcr, dat.fit.bcr, xlim = c(0, 0.05))
vari_plot_10

#fitting the kriged data with variogram and detrended data
f.bcr <- as.formula(cd_10 ~ X + Y) 

dat.krg <- krige(f.bcr, bcr, grd, dat.fit.bcr)

r <- raster(dat.krg)
r.m <- mask(r, base_spatial)

#plot!
i_10 <- 
tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Crowd Density at 10pm", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
  tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants'))+
  tm_legend(legend.outside=TRUE) 
i_10

#11pm

f.bcr <- as.formula(cd_11 ~ lon + lat) 

var.smpl.bcr <- variogram(f.bcr, bcr, cloud = FALSE, cutoff=1, width=0.0025)

dat.fit.bcr  <- fit.variogram(var.smpl.bcr, fit.ranges = F, fit.sills = F,
                          vgm(psill=350, model="Wav", range=0.025, nugget=800))

# The following plot allows us to assess the fit
vari_plot_11 <- plot(var.smpl.bcr, dat.fit.bcr, xlim = c(0, 0.05))
vari_plot_11


f.bcr <- as.formula(cd_11 ~ X + Y) 

dat.krg <- krige(f.bcr, bcr, grd, dat.fit.bcr)

r <- raster(dat.krg)
r.m <- mask(r, base_spatial)

#tm_shape(base_map) + tm_borders(col='black') + tm_fill(col = "white") +
i_11 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Crowd Density at 11pm", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
    tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants'))+
  tm_legend(legend.outside=TRUE) 
i_11

#12am

f.bcr <- as.formula(cd_12 ~ lon + lat) 

var.smpl.bcr <- variogram(f.bcr, bcr, cloud = FALSE, cutoff=1, width=0.0025)

dat.fit.bcr  <- fit.variogram(var.smpl.bcr, fit.ranges = F, fit.sills = F,
                          vgm(psill=250, model="Wav", range=0.025, nugget=900))

# The following plot allows us to assess the fit
vari_plot_12 <- plot(var.smpl.bcr, dat.fit.bcr, xlim = c(0, 0.05))
vari_plot_12


f.bcr <- as.formula(cd_12 ~ X + Y) 

dat.krg <- krige(f.bcr, bcr, grd, dat.fit.bcr)

r <- raster(dat.krg)
r.m <- mask(r, base_spatial)

#tm_shape(base_map) + tm_borders(col='black') + tm_fill(col = "white") +
i_12 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Crowd Density at 12am", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
    tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants'))+
  tm_legend(legend.outside=TRUE) 
i_12

#1am

f.bcr <- as.formula(cd_1 ~ lon + lat) 

var.smpl.bcr <- variogram(f.bcr, bcr, cloud = FALSE, cutoff=1, width=0.0025)

dat.fit.bcr  <- fit.variogram(var.smpl.bcr, fit.ranges = F, fit.sills = F,
                          vgm(psill=50, model="Wav", range=0.01, nugget=900))

# The following plot allows us to assess the fit
vari_plot_1 <- plot(var.smpl.bcr, dat.fit.bcr, xlim = c(0, 0.05))
vari_plot_1


f.bcr <- as.formula(cd_1 ~ X + Y) 

dat.krg <- krige(f.bcr, bcr, grd, dat.fit.bcr)

r <- raster(dat.krg)
r.m <- mask(r, base_spatial)

#tm_shape(base_map) + tm_borders(col='black') + tm_fill(col = "white") +
i_1 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Crowd Density at 1am", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
    tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants'))+
  tm_legend(legend.outside=TRUE) 
i_1

#2am

f.bcr <- as.formula(cd_2 ~ lon + lat) 

var.smpl.bcr <- variogram(f.bcr, bcr, cloud = FALSE, cutoff=1, width=0.0025)

dat.fit.bcr  <- fit.variogram(var.smpl.bcr, fit.ranges = F, fit.sills = F,
                          vgm(psill=150, model="Wav", range=0.015, nugget=550))

# The following plot allows us to assess the fit
vari_plot_2 <- plot(var.smpl.bcr, dat.fit.bcr, xlim = c(0, 0.05))
vari_plot_2


f.bcr <- as.formula(cd_2 ~ X + Y) 

dat.krg <- krige(f.bcr, bcr, grd, dat.fit.bcr)

r <- raster(dat.krg)
r.m <- mask(r, base_spatial)

#tm_shape(base_map) + tm_borders(col='black') + tm_fill(col = "white") +
i_2 <- tm_shape(r.m) + 
  tm_raster(n=13, stretch.palette=FALSE, 
            title="Crowd Density at 2am", palette = 'Reds', style = 'fixed',breaks = c(0,10,20,30,40,50,60,70,80,90,100,110,120)) +
    tm_shape(base_map) + tm_borders(col='black') + tm_fill(alpha = 0) +
  tm_shape(bcr) + tm_dots(size=0.1, col = 'cat', palette = c('green2','blue'), title = 'Label', label = c('Nightlife locations','Late night restaurants'))+
  tm_legend(legend.outside=TRUE) 
i_2


```


### Creating Animation!

```{r spatial animations}


tmap_save(i_10, filename = "plots/crowd_density/spatial trend at 10pm.png")
tmap_save(i_11, filename = "plots/crowd_density/spatial trend at 11pm.png")
tmap_save(i_12, filename = "plots/crowd_density/spatial trend at 12am.png")
tmap_save(i_1, filename = "plots/crowd_density/spatial trend at 1am.png")
tmap_save(i_2, filename = "plots/crowd_density/spatial trend at 2am.png")

img_list <- image_read(c("plots/crowd_density/spatial trend at 10pm.png","plots/crowd_density/spatial trend at 11pm.png","plots/crowd_density/spatial trend at 12am.png","plots/crowd_density/spatial trend at 1am.png","plots/crowd_density/spatial trend at 2am.png"))
animation <- image_animate(img_list, fps = 0.5)  # Adjust fps as needed
image_write(animation, path = "plots/crowd_density/spatial trends.gif")


```



## 5. KDE to recommend new late night restaurants

```{r kde-of-nightlife}
base_map <- st_transform(central_area, crs=st_crs(central_area)) #Reset crs

nightlife_points <- as(nightlife_sf, "Spatial")
nightlife_centers <- kde.points(nightlife_points, h = 500)
nightlife_centers_sf <- st_as_sf(nightlife_centers)

tmap_mode("plot")
tm_shape(nightlife_centers) + tm_raster()
```


```{r get-kde-polygons}
reclass_values <- c(0,0.0000002,1, #reclassify kde values from 0-0.0000002 in group 1 and so on
                    0.0000002,0.0000004,2,
                    0.0000004,0.0000006,3,
                    0.0000006,0.0000008,4,
                    0.0000008,0.0000010,5,
                    0.0000010,0.0000012,6)

reclass_nightlife_centers <- reclassify(as(nightlife_centers, "RasterLayer"), reclass_values) #reclassify kde values to groups
nightlife_centers_poly <- rasterToPolygons(reclass_nightlife_centers, dissolve = T) #to make a polygon layer
nightlife_centers_poly <- st_as_sf(nightlife_centers_poly) #to make an SF object
nightlife_centers_poly <- nightlife_centers_poly[-c(1),] #remove polys with low kde values 
nightlife_centers_poly <- st_cast(nightlife_centers_poly,'POLYGON') #to split multipolygon to polygon to obtain centers
```

```{r plot-kde}
nightlife_centers_poly <- st_transform(nightlife_centers_poly, crs=st_crs(subzones)) # project to basemap's CRS

# Plot of nightlife location and density polygons
tm_shape(base_map) + tm_borders() + 
  tm_shape(nightlife_centers_poly) + tm_fill(col = 'kde') + 
  tm_shape(nightlife_sf) + tm_bubbles(size = 0.001, col = "black")  +
  tm_compass(type="8star", size = 2) +
  tm_scale_bar(width = 0.15) +
  tm_layout(legend.format = list(digits = 0),
            legend.position = c("left", "top"),
            legend.text.size = 0.5, 
            legend.title.size = 1,
            title="KDE of Nightlife locations in Central Singapore",
            title.position = c('left', 'top'))
```

```{r get restaurants-in-kde-region}
# Filter for nightlife locations within KDE regions
all_kde <- st_union(nightlife_centers_poly)
nightlife_kde <- st_intersection(all_kde, nightlife_sf)
print(length(nightlife_kde))
print(nrow(nightlife_sf))
print(sprintf("A total of %d nightlife locations lie in the filtered KDE regions, out of %d original nightlife locations. That means %f of nightlife locations lie in the filtered KDE regions.", length(nightlife_kde), nrow(nightlife_sf), length(nightlife_kde)/nrow(nightlife_sf)))
```
```{r recommended-restaurants}

midnight_restaurants <- restaurants_sf %>%
  filter(midnight_restaurant == "True")

non_midnight_restaurants <- restaurants_sf %>%
  filter(midnight_restaurant == "False")

# Combine polygons by KDE class and find its intersection with restaurants
grouped_union <- nightlife_centers_poly %>%
  group_by(kde) %>%
  summarise(geometry = st_union(geometry))

nightlife_intersection <- st_intersection(grouped_union, nightlife_sf)

kde_nightlife_count <- nightlife_intersection %>%
  group_by(kde) %>%
  summarise(nightlife_count = n())

print(kde_nightlife_count)

# Combine polygons by KDE class and find its intersection with restaurants
restaurant_intersection <- st_intersection(grouped_union, non_midnight_restaurants)
midnight_restaurant_intersection <- st_intersection(grouped_union, midnight_restaurants)

kde_restaurant_count <- restaurant_intersection %>%
  group_by(kde) %>%
  summarise(restaurant_count = n())

# View the result
print(kde_restaurant_count)
```

```{r}
# Plot of all non-midnight restaurants within the extracted polygons
tm_shape(restaurant_intersection) +
  tm_bubbles(size = 0.1, col = "kde") +
  tm_shape(midnight_restaurant_intersection) +
  tm_bubbles(size = 0.1, col = "black") +
  tm_compass(type="8star", size = 2) +
  tm_scale_bar(width = 0.15) +
  tm_layout(legend.format = list(digits = 0),
            legend.position = c("left", "top"),
            legend.text.size = 0.5, 
            legend.title.size = 1,
            title="Recommended Restaurants in Central Singapore",
            title.position = c('left', 'top'))
```