example_sfnetworks.md

# packages
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sfnetworks)
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter

# data
download.file("https://github.com/JovaniSouza/JovaniSouza5/raw/master/Test.zip", "Test.zip")
download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
unzip("Test.zip")

roads = st_read("Test/regionbrazil.shp", quiet = TRUE)
points = st_read("Example/Points/Points.shp", quiet = TRUE)
roads_trf = st_transform(roads, st_crs(points)) %>% 
  st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant

# build sfnetwork
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

# Estimate shortest path between 3rd point and 4th point
st_shortest_paths(net, points[3, ], points[4, ])
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> $vpath
#> $vpath[[1]]
#> + 1/98316 vertex, from 1f5f38c:
#> [1] 65808
#> 
#> 
#> $epath
#> NULL
#> 
#> $predecessors
#> NULL
#> 
#> $inbound_edges
#> NULL

# It looks like there is not path, but why? The problem is that both points are
# associated to the same node:
st_nearest_feature(points[3, ], net %>% activate(nodes) %>% st_as_sf())
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> [1] 65808
st_nearest_feature(points[4, ], net %>% activate(nodes) %>% st_as_sf())
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> [1] 65808

# So it's not a problem of the function of the package but simply both points
# are matched with the same node so there is no "shortest path".

^{Created on 2020-07-12 by the reprex package (v0.3.0)}

Thanks @agila5. I updated the packages mentioned, and it decreased the time considerably. Thanks for the tip.

I will wait for the example. If it is not difficult, you could use these same roads and points of the code above.

Thanks again!

Here it is:

# packages
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sfnetworks)
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(tmap)
library(tictoc)

# data
download.file("https://github.com/JovaniSouza/JovaniSouza5/raw/master/Test.zip", "Test.zip")
download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
unzip("Test.zip")

roads = st_read("Test/regionbrazil.shp", quiet = TRUE)
points = st_read("Example/Points/Points.shp", quiet = TRUE)
roads_trf = st_transform(roads, st_crs(points)) %>% 
  st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant

# build sfnetwork
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

# routing
from = c(-49.95058, -24.77502)
to = c(-49.91084, -24.75200)
p1 = sf::st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = sf::st_crs(net))
p2 = sf::st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = sf::st_crs(net))

r = tidygraph::convert(net, to_spatial_shortest_paths, p1, p2)
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
r2 = tidygraph::convert(net, to_spatial_shortest_paths, points[30, ], points[8, ])
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar

# Extract the bbox for r and r2
bbox_r = st_as_sfc(r %>% activate(edges) %>% st_bbox())
bbox_r2 = st_as_sfc(r2 %>% activate(edges) %>% st_bbox())

# merge the two bboxes and make them a little bit bigger
merged_bboxes = st_as_sfc(st_bbox(st_union(bbox_r, bbox_r2)))
#> although coordinates are longitude/latitude, st_union assumes that they are planar

# filter the net
small_net = st_filter(net, merged_bboxes)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

# plot
tic()
tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
    tm_lines() + 
tm_shape(rbind(p1, p2)) + 
  tm_dots(col = "red", size = 0.5) + 
tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "red", lwd = 3) + 
tm_shape(points[c(30, 8), ]) + 
  tm_dots(col = "blue", size = 0.5) + 
tm_shape(r2 %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "blue", lwd = 3)

toc()
#> 0.32 sec elapsed

# We can also make a little bit bigger bbox (note that you are working with
# long/lat coordinates so it's no easy to specify an expanding factor)
merged_bboxes = st_as_sfc(st_bbox(st_union(bbox_r, bbox_r2)) * c(1.0005, 1.002, 0.998, 0.998))
#> although coordinates are longitude/latitude, st_union assumes that they are planar

# filter the net
small_net = st_filter(net, merged_bboxes)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

# plot
tic()
tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
  tm_lines() + 
  tm_shape(rbind(p1, p2)) + 
  tm_dots(col = "red", size = 0.5) + 
  tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "red", lwd = 3) + 
  tm_shape(points[c(30, 8), ]) + 
  tm_dots(col = "blue", size = 0.5) + 
  tm_shape(r2 %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "blue", lwd = 3)

toc()
#> 0.33 sec elapsed

^{Created on 2020-07-14 by the reprex package (v0.3.0)}

Very good @agila5!

It's much faster than before. Even using the complete_roads shapefile, that file that has multiple roads. It helped me a lot, I appreciate your help.

I believe that I have no more questions.

I am happy to help too, should you need it.

Best Regards.

Hello @agila5. Just a quick question: If I didn't want to use my points file (that is in shapefile) but a data frame. I want to see the distance between property 1 and 7 from my df database, but without directly using the coordinates in from and to. Is it possible to do it in some other way? For example use for from = c(df[1,c("Longitude")], df[1,c("Latitude")]) and to = c(df[7,c("Longitude")], df[7,c("Latitude")]).

library(sf)
library(sfnetworks)
library(tidygraph)
library(tmap)

  df<-structure(list(Property = c(1,2,3,4,5,6,7), 
                     Latitude = c(-24.779225, -24.789635, -24.763461, -24.794394, -24.747102,-24.781307,-24.761081), 
                     Longitude = c(-49.934816, -49.922324, -49.911616, -49.906262, -49.890796,-49.8875254,-49.8875254), 
                     Waste = c(526, 350, 526, 469, 285, 433, 456)), class = "data.frame", row.names = c(NA, -7L))

roads = st_read("Test/regionbrazil.shp", quiet = TRUE)
points = st_read("Example/Points/Points.shp", quiet = TRUE)
roads_trf = st_transform(roads, st_crs(points)) %>% 
  st_cast("LINESTRING")

# build sfnetwork
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

# routing
from = c(-49.87, -24.80)
to = c(-49.92, -24.78)
p1 = sf::st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = sf::st_crs(net))
p2 = sf::st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = sf::st_crs(net))
r = tidygraph::convert(net, to_spatial_shortest_paths, p1, p2)

# Extract the bbox for r 
bbox_r = st_as_sfc(r %>% activate(edges) %>% st_bbox())

# filter the net
small_net = st_filter(net, bbox_r)

# plot
tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
  tm_lines() + 
  tm_shape(rbind(p1, p2)) + 
  tm_dots(col = "red", size = 0.5) + 
  tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "red", lwd = 3)

Is it possible to do it in some other way? For example use for from = c(df[1,c("Longitude")], df[1,c("Latitude")]) and to = c(df[7,c("Longitude")], df[7,c("Latitude")])

Yes, it should work. Did you test it?

Hi, @agila5. Thanks for the reply.

This part I did works, but how does that part look?

roads_trf = st_transform (roads, st_crs (points))%>% st_cast ("LINESTRING")

I don't want to use the points file, just the df database and the roads file. So how would it look?

Best regards.

The function st_transform() is used to change the CRS of roads object and you cannot change it CRS to st_crs(points) unless you load the points data. Can you really skip that step?

But then I would need to work on that part too, right? because it calling roads_trf.

net = as_sfnetwork(roads_trf, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

I replaced only with roads, but it doesn't work, it gives error. Please could you test it on yours? see below how I did it:

library(sf)
library(sfnetworks)
library(tidygraph)
library(tmap)

download.file("https://github.com/JovaniSouza/JovaniSouza5/raw/master/Test.zip", "Test.zip")
unzip("Test.zip")

df<-structure(list(Property = c(1,2,3,4,5,6,7), 
                   Latitude = c(-24.779225, -24.789635, -24.763461, -24.794394, -24.747102,-24.781307,-24.761081), 
                   Longitude = c(-49.934816, -49.922324, -49.911616, -49.906262, -49.890796,-49.8875254,-49.8875254), 
                   Waste = c(526, 350, 526, 469, 285, 433, 456)), class = "data.frame", row.names = c(NA, -7L))

roads = st_read("Test/regionbrazil.shp", quiet = TRUE)

# build sfnetwork
net = as_sfnetwork(roads, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

# routing
from = c(df[1,c("Longitude")], df[1,c("Latitude")])
to = c(df[7,c("Longitude")], df[7,c("Latitude")])
p1 = sf::st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = sf::st_crs(net))
p2 = sf::st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = sf::st_crs(net))
r = tidygraph::convert(net, to_spatial_shortest_paths, p1, p2)

# Extract the bbox for r 
bbox_r = st_as_sfc(r %>% activate(edges) %>% st_bbox())

# filter the net
small_net = st_filter(net, bbox_r)

# plot
tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
  tm_lines() + 
  tm_shape(rbind(p1, p2)) + 
  tm_dots(col = "red", size = 0.5) + 
  tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "red", lwd = 3)

Hi. The error is

Error in as_sfnetwork.sf(roads, directed = FALSE) : 
  Only geometries of type LINESTRING or POINT are allowed

and, IMO, it's pretty clear. You need to transform the geometry of the roads object into LINESTRING(s). The following is a working example:

# packages
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sfnetworks)
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(tmap)

# data
download.file("https://github.com/JovaniSouza/JovaniSouza5/raw/master/Test.zip", "Test.zip")
unzip("Test.zip")

df <- structure(
  list(
    Property = c(1,2,3,4,5,6,7), 
    Latitude = c(-24.779225, -24.789635, -24.763461, -24.794394, -24.747102,-24.781307,-24.761081),
    Longitude = c(-49.934816, -49.922324, -49.911616, -49.906262, -49.890796,-49.8875254,-49.8875254), 
    Waste = c(526, 350, 526, 469, 285, 433, 456)
  ), 
  class = "data.frame", 
  row.names = c(NA, -7L)
)

roads = st_read("Test/regionbrazil.shp", quiet = TRUE) %>% 
  st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant

# build sfnetwork
net = as_sfnetwork(roads, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

# routing
from = c(df[1, c("Longitude")], df[1, c("Latitude")])
to = c(df[7, c("Longitude")], df[7, c("Latitude")])
p1 = st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = st_crs(net))
p2 = st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = st_crs(net))
r = tidygraph::convert(net, to_spatial_shortest_paths, p1, p2)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar

# Extract the bbox for r 
bbox_r = st_as_sfc(r %>% activate(edges) %>% st_bbox())

# filter the net
small_net = st_filter(net, bbox_r)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

# plot
tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
  tm_lines() + 
  tm_shape(rbind(p1, p2)) + 
  tm_dots(col = "red", size = 0.5) + 
  tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "red", lwd = 3)

^{Created on 2020-07-18 by the reprex package (v0.3.0)}

You need to be 100% sure that the CRS of the roads object is the same CRS as the points otherwise you can get misleading results. Anyway, I think you should read the vignettes of sf and the introductory vignette of sfnetworks and maybe the first 7 chapters of this book to understand a bit more details on how to manipulate sf objects.

Thanks again @agila5 for the help. Sorry if I asked a few simple questions, is that this area of programming in R is new to me, to be more precise I have been working in R for just three months. Thanks for the tips pointed out and I will definitely read your recommendations above.
Best Regards.

No worries and good luck with your project. If you have any question on sfnetworks feel free to ask.

Thank you very much @agila5.

Hello @agila5, I hope you are well. I have a question regarding my project that involves the sfnetoworks package and Shiny format, I don't know your knowledge about shiny. But I will insert below, if you have any ideas or solutions, I would appreciate it.

I am inserting three executable codes below, the first generates a map using sftnetworks package, showing the route between two locations. In this case, the two locations to generate the map were defined: from = c(df_spec_clust[1, c("Longitude")], df_spec_clust[1, c("Latitude")])
and to = c (df_spec_prop [4, c ("Longitude")], df_spec_prop [4, c ("Latitude")])]. In the second, I would like to generate the map in Shiny format, but without defining the locations exactly as I did in the first code. I would like them to be selected from the filters I created (Filter 1 and Filter 2). However, I am unable to generate the map. Could you help me ? To show you I managed to generate the map correctly in the third code for the problem in question but using another package (leaflet). However, I still couldn't think of a way to make it work using the sfnetworks package. Any help is appreciated.

Thank you!

First code

library(sf)
library(sfnetworks)
library(tmap)
library(rdist)
library(geosphere)

#for the roads file
download.file("https://github.com/JovaniSouza/JovaniSouza5/raw/master/Test.zip", "Test.zip")
unzip("Test.zip")

#database df
df <- structure(
  list(Property = c(1,2,3,4,5,6,7), Latitude = c(-24.779225, -24.789635, -24.763461, -24.794394, -24.747102,-24.781307,-24.761081),
    Longitude = c(-49.934816, -49.922324, -49.911616, -49.906262, -49.890796,-49.8875254,-49.8875254), 
    Waste = c(526, 350, 526, 469, 285, 433, 456)),class = "data.frame", row.names = c(NA, -7L))

#clusters
coordinates<-df[c("Latitude","Longitude")]
d<-as.dist(distm(coordinates[,2:1]))
fit.average<-hclust(d,method="average") 
k=3
clusters<-cutree(fit.average, k) 
nclusters<-matrix(table(clusters))  
df$cluster <- clusters 

#Create database df1
center<-matrix(nrow=k,ncol=2)
for(i in 1:k){
  center[i,]<-c(weighted.mean(subset(df,cluster==i)$Latitude,subset(df,cluster==i)$Waste),
                     weighted.mean(subset(df,cluster==i)$Longitude,subset(df,cluster==i)$Waste))}
coordinates$cluster<-clusters 
center<-cbind(center,matrix(c(1:k),ncol=1)) 
df1<-as.data.frame(center)
colnames(df1) <-c("Latitude", "Longitude", "cluster")

#specific cluster and specific property
df_spec_clust <- df1[df1$cluster,]
df_spec_prop<-df[df$Property,]

#create map
roads = st_read("Test/regionbrazil.shp", quiet = TRUE) %>% 
  st_cast("LINESTRING")

# build sfnetwork
net = as_sfnetwork(roads, directed = FALSE) %>%
  activate("edges") %>%
  dplyr::mutate(weight = edge_length())

# routing
from = c(df_spec_clust[1, c("Longitude")], df_spec_clust[1, c("Latitude")])
to = c(df_spec_prop[4, c("Longitude")], df_spec_prop[4, c("Latitude")])
p1 = st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = st_crs(net))
p2 = st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = st_crs(net))
r = tidygraph::convert(net, to_spatial_shortest_paths, p1, p2)

# Extract the bbox for r 
bbox_r = st_as_sfc(r %>% activate(edges) %>% st_bbox())


# filter the net
small_net = st_filter(net, bbox_r)

# plot
plot1<-tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
  tm_lines() + 
  tm_shape(rbind(p1, p2)) + 
  tm_dots(col = "red", size = 0.5) + 
  tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
  tm_lines(col = "red", lwd = 3)
plot1

Map generated by the code above

Second code

library(shiny)
library(rdist)
library(geosphere)
library(shinythemes)
library(sf)
library(tidygraph)
library(sfnetworks)
library(tmap)

#for the roads file
 download.file("https://github.com/JovaniSouza/JovaniSouza5/raw/master/Test.zip", "Test.zip")
 unzip("Test.zip")

function.cl<-function(df,k,Filter1,Filter2){
  
  #database df
  df <- structure(
    list(Property = c(1,2,3,4,5,6,7), Latitude = c(-24.779225, -24.789635, -24.763461, -24.794394, -24.747102,-24.781307,-24.761081),
         Longitude = c(-49.934816, -49.922324, -49.911616, -49.906262, -49.890796,-49.8875254,-49.8875254), 
         Waste = c(526, 350, 526, 469, 285, 433, 456)),class = "data.frame", row.names = c(NA, -7L))
  
  #clusters
  coordinates<-df[c("Latitude","Longitude")]
  d<-as.dist(distm(coordinates[,2:1]))
  fit.average<-hclust(d,method="average") 
  clusters<-cutree(fit.average, k) 
  nclusters<-matrix(table(clusters))  
  df$cluster <- clusters 
  
  #Create database df1
  center<-matrix(nrow=k,ncol=2)
  for(i in 1:k){
    center[i,]<-c(weighted.mean(subset(df,cluster==i)$Latitude,subset(df,cluster==i)$Waste),
                  weighted.mean(subset(df,cluster==i)$Longitude,subset(df,cluster==i)$Waste))}
  coordinates$cluster<-clusters 
  center<-cbind(center,matrix(c(1:k),ncol=1)) 
  df1<-as.data.frame(center)
  colnames(df1) <-c("Latitude", "Longitude", "cluster")
 
  # specific cluster and specific property
  df_spec_clust <- df1[df1$cluster==Filter1,]
  df_spec_prop<-df[df$Property==Filter2,]
  
 
  #create map
 
  roads = st_read("Test/regionbrazil.shp", quiet = TRUE) %>% 
    st_cast("LINESTRING")
  
  # build sfnetwork
  net = as_sfnetwork(roads, directed = FALSE) %>%
    activate("edges") %>%
    dplyr::mutate(weight = edge_length())
  
  # routing
  from = c(df_spec_clust[1, c("Longitude")], df_spec_clust[1, c("Latitude")])
  to = c(df_spec_prop[4, c("Longitude")], df_spec_prop[4, c("Latitude")])
  p1 = st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = st_crs(net))
  p2 = st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = st_crs(net))
  r = tidygraph::convert(net, to_spatial_shortest_paths, p1, p2)
  
  # Extract the bbox for r 
  bbox_r = st_as_sfc(r %>% activate(edges) %>% st_bbox())
  
  
  # filter the net
  small_net = st_filter(net, bbox_r)
  
  # plot
  plot1<-tm_shape(small_net %>% activate(edges) %>% st_as_sf()) +
    tm_lines() + 
    tm_shape(rbind(p1, p2)) + 
    tm_dots(col = "red", size = 0.5) + 
    tm_shape(r %>% activate(edges) %>% st_as_sf()) + 
    tm_lines(col = "red", lwd = 3)

  return(list(
    "Plot1" = plot1,
    "Data" =  df
  ))
}

ui <- bootstrapPage(
  navbarPage(theme = shinytheme("flatly"), collapsible = TRUE,
             "Cl", 
          tabPanel("",
           sidebarLayout(
             sidebarPanel(
               sliderInput("Slider", h5(""),
                           min = 2, max = 4, value = 3),
               selectInput("Filter1", label = h4("Select just one cluster"),""),
               selectInput("Filter2",label=h4("Select the cluster property"),""),
             ),
             mainPanel(
               tabsetPanel(
                 tabPanel("Map", plotOutput("Map1"))))
           ))))

server <- function(input, output, session) {
  
  Modelcl<-reactive({
    function.cl(df,input$Slider,input$Filter1,input$Filter2)
  })
  

  output$Map1 <- renderPlot({
    Modelcl()[[1]]
  })
  
  observeEvent(input$Slider, {
    abc <- req(Modelcl()$Data)
    updateSelectInput(session,'Filter1',
                      choices=sort(unique(abc$cluster)))
  }) 
  
  observeEvent(input$Filter1,{
    abc <- req(Modelcl()$Data) %>% filter(cluster == as.numeric(input$Filter1))
    updateSelectInput(session,'Filter2',
                      choices=sort(unique(abc$Property)))
  }) 
  
  
}

shinyApp(ui = ui, server = server)

Map generated but using leaflet package (It works)

library(shiny)
library(rdist)
library(geosphere)
library(shinythemes)
library(leaflet)
library(tidygraph)

function.cl<-function(df,k,Filter1,Filter2){
  
  #database df
  df <- structure(
    list(Property = c(1,2,3,4,5,6,7), Latitude = c(-24.779225, -24.789635, -24.763461, -24.794394, -24.747102,-24.781307,-24.761081),
         Longitude = c(-49.934816, -49.922324, -49.911616, -49.906262, -49.890796,-49.8875254,-49.8875254), 
         Waste = c(526, 350, 526, 469, 285, 433, 456)),class = "data.frame", row.names = c(NA, -7L))
  
  #clusters
  coordinates<-df[c("Latitude","Longitude")]
  d<-as.dist(distm(coordinates[,2:1]))
  fit.average<-hclust(d,method="average") 
  clusters<-cutree(fit.average, k) 
  nclusters<-matrix(table(clusters))  
  df$cluster <- clusters 
  
  #Create database df1
  center<-matrix(nrow=k,ncol=2)
  for(i in 1:k){
    center[i,]<-c(weighted.mean(subset(df,cluster==i)$Latitude,subset(df,cluster==i)$Waste),
                  weighted.mean(subset(df,cluster==i)$Longitude,subset(df,cluster==i)$Waste))}
  coordinates$cluster<-clusters 
  center<-cbind(center,matrix(c(1:k),ncol=1)) 
  df1<-as.data.frame(center)
  colnames(df1) <-c("Latitude", "Longitude", "cluster")
  
  #specify cluster and specific cluster and specific propertie
  df_spec_clust <- df1[df1$cluster==Filter1,]
  df_spec_prop<-df[df$Property==Filter2,]
  
  
  #color for map
  ai_colors <-c("red","gray","blue","orange","green","beige","darkgreen","lightgreen", "lightred", "darkblue","lightblue",
                "purple","darkpurple","pink", "cadetblue","white","darkred", "lightgray","black")
  clust_colors <- ai_colors[df$cluster]
  icons <- awesomeIcons(
    icon = 'ios-close',
    iconColor = 'black',
    library = 'ion',
    markerColor =  clust_colors)
  
  # create icon for map
  leafIcons <- icons(
    iconUrl = ifelse(df1$cluster,
                     
                     "https://image.flaticon.com/icons/svg/542/542461.svg"
    ),
    iconWidth = 30, iconHeight = 40,
    iconAnchorX = 25, iconAnchorY = 12)
  
  html_legend <- "<img src='https://image.flaticon.com/icons/svg/542/542461.svg'>"
  
# create map
  if(nrow(df_spec_clust)>0){
    clust_colors <- ai_colors[df_spec_clust$cluster]
    icons <- awesomeIcons(
      icon = 'ios-close',
      iconColor = 'black',
      library = 'ion',
      markerColor =  clust_colors)
    
  m1<-leaflet(df_spec_clust) %>% addTiles() %>% 
    addMarkers(~Longitude, ~Latitude, icon = leafIcons) %>%
    addAwesomeMarkers(leaflet(df_spec_prop) %>% addTiles(), lat=~df_spec_prop$Latitude, lng = ~df_spec_prop$Longitude, icon= icons,label=~cluster)

  for(i in 1:nrow(df_spec_clust)){
    df_line <- rbind(df_spec_prop[,c("Latitude","Longitude")],
                     df_spec_clust[i,c("Latitude","Longitude")])
    m1 <- m1 %>%
      addPolylines(data = df_line,
                   lat=~Latitude,
                   lng = ~Longitude,
                   color="red")
  }
  plot1<-m1} else plot1 <- NULL

  return(list(
    "Plot1" = plot1,
    "Data"= df
  ))
}

ui <- bootstrapPage(
  navbarPage(theme = shinytheme("flatly"), collapsible = TRUE,
             "Cl", 
             tabPanel("",
                      sidebarLayout(
                        sidebarPanel(
                          sliderInput("Slider", h5(""),
                                      min = 2, max = 4, value = 3),
                          selectInput("Filter1", label = h4("Select just one cluster"),""),
                          selectInput("Filter2",label=h4("Select the cluster property"),""),
                        ),
                        mainPanel(
                          tabsetPanel(
                            tabPanel("Map", uiOutput("Map1"))))
                      ))))

server <- function(input, output, session) {
  
  Modelcl<-reactive({
    function.cl(df,input$Slider,input$Filter1,input$Filter2)
  })
  
  output$Map1 <- renderUI({ 
    if(input$Filter1!="") 
      leafletOutput("Leaf1",width = "95%", height = "600") })

  output$Leaf1 <- renderLeaflet({
    req(Modelcl())[[1]]
  })
  
  
  observeEvent(input$Slider, {
    abc <- req(Modelcl()$Data)
    updateSelectInput(session,'Filter1',
                      choices=sort(unique(abc$cluster)))
  }) 
  
  observeEvent(input$Filter1,{
    abc <- req(Modelcl()$Data) %>% filter(cluster == as.numeric(input$Filter1))
    updateSelectInput(session,'Filter2',
                      choices=sort(unique(abc$Property)))
  }) 
  
  
}

shinyApp(ui = ui, server = server)

Hi @JovaniSouza, I'm sorry but I've never used shiny apart from the trivial examples so I cannot help you. I can still check the code in the next days but I'm not sure I will be able to solve the problem.

Probably it's easier if you can create a smaller example showing the problem.

Thanks for the reply @agila5. Do you find it interesting to include this issue in the sfnetworks package issues?

I don't think it's a problem with sfnetworks since the non-shiny-app approach works. IMO the best way for you would be creating a smaller example and add your question to StackOverflow (or similar websites)

Ah understood. Anyway, if you can think of any possibility that I can resolve the above question, please do not hesitate to contact me.

Thanks again @agila5.