Test Agh Food
Matt Broerman
10/18/2021
Background
This demo is adapted from Edgar Ruiz’s RStudio conference.
Facilities located within Allegheny County that produce, distribute and sell food products are subject to mandatory, routine inspection by one of the health department’s staff of environmental health specialists.
Azure SQL Server
con <-
dbConnect(
odbc::odbc(),
Driver = "ODBC Driver 17 for SQL Server",
Server = Sys.getenv("paccmdb_server"),
Database = Sys.getenv("paccmdb_pulce"),
Uid = Sys.getenv("pitt_userid"),
Pwd = Sys.getenv("pitt_pwd"),
Encrypt = "yes",
TrustServerCertificate = "no",
Authentication = "ActiveDirectoryPassword",
timeout = 10
)dplyr with SQL Server
Instead of importing the entire data set, we just create a pointer to the table inside SQL Server:
service <- tbl(con, "test_agh_food")Standard dplyr commands work inside SQL Server:
glimpse(service)## Rows: ??
## Columns: 23
## Database: Microsoft SQL Server 12.00.2195[MJB357@pitt.edu@paccmdb/paccmdb]
## $ id <dbl> 10762, 10781, 10693, 10711, 10734, 10655, 10664, 10668, …
## $ facility_name <chr> "North Versailles VFD / South Wilmerding Social Club", "…
## $ num <chr> "830", "128-30", "3113", "540", "415", "5775", "2425", "…
## $ street <chr> "Sylvan Ave", "Ormsby Ave", "W Liberty Ave", "3rd Avenu…
## $ city <chr> "North Versailles", "Pittsburgh", "Pittsburgh", "Bridgev…
## $ state <chr> "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "PA", "P…
## $ zip <dbl> 15137, 15210, 15216, 15017, 15239, 15102, 15102, 15017, …
## $ municipal <chr> "North Versailles", "Mt Oliver", "Dormont", "Bridgeville…
## $ category_cd <dbl> 250, 250, 201, 201, 408, 212, 201, 211, 201, 201, 250, 2…
## $ description <chr> "Social Club-Bar Only", "Social Club-Bar Only", "Restaur…
## $ p_code <dbl> 1, 7, 2, 1, 1, 7, 1, 7, 1, 1, 2, 1, 7, 0, 7, 2, 1, 1, 2,…
## $ fdo <date> 1969-01-01, 1984-06-17, 1984-06-17, 1984-06-17, 1957-01…
## $ bus_st_date <date> 1969-01-01, 1984-06-17, 1984-06-17, 1984-06-17, 1957-01…
## $ bus_cl_date <date> NA, NA, NA, NA, NA, NA, 2017-09-28, NA, NA, NA, 2017-07…
## $ seat_count <dbl> 60, 35, 41, 75, 125, 13, 70, 0, 30, 30, 30, 220, 30, 45,…
## $ noroom <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ sq_feet <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ status <chr> "Active", "Out of Business", "Out of Business", "Active"…
## $ placard_st <dbl> 1, NA, 1, 1, 1, NA, 1, NA, 1, 1, 1, 1, NA, 1, NA, NA, 1,…
## $ x <dbl> -79.79541, -79.86383, -80.03935, -80.13519, -79.70589, -…
## $ y <dbl> 40.38935, 40.42811, 40.39097, 40.35593, 40.46730, 40.348…
## $ address <chr> "830 Sylvan Ave North Versailles, PA 15137", "128-30 Orm…
## $ year <dbl> 1969, 1984, 1984, 1984, 1957, 1999, 1989, 1984, 1984, 19…tally(service)## # Source: lazy query [?? x 1]
## # Database: Microsoft SQL Server 12.00.2195[MJB357@pitt.edu@paccmdb/paccmdb]
## n
## <int>
## 1 24357Which are the most common types?
service %>%
count(description, sort = TRUE)## # Source: lazy query [?? x 2]
## # Database: Microsoft SQL Server 12.00.2195[MJB357@pitt.edu@paccmdb/paccmdb]
## # Ordered by: desc(n)
## description n
## <chr> <int>
## 1 Restaurant without Liquor 4978
## 2 Restaurant with Liquor 3624
## 3 Chain Restaurant without Liquor 3082
## 4 Retail/Convenience Store 1262
## 5 Chain Retail/Convenience Store 1114
## 6 Chain Packaged Food Only 1027
## 7 Packaged Food Only 948
## 8 Mobile â?? Tier I (Unopened Pre-packaged Food Only) 929
## 9 Chain Restaurant with Liquor 822
## 10 Child Food Service 501
## # … with more rowsSome data cleaning…
# service <- service %>%
# filter(latitude > 10)Visualizations using dbplot, ggplot2 and leaflet
Easily view the location of the facilities without importing all of the data.
library(dbplot)
restaurant <- "Restaurant without Liquor"
service %>%
filter(description == restaurant) %>%
dbplot_raster(x, y)
Use the new db_computer_raster2() function to obtain the pre-calculated squares, with their limits included.
restaurants <- service %>%
filter(description == restaurant) %>%
db_compute_raster2(x, y, resolution = 30) %>%
mutate(of_max = `n()` / max(`n()`))
restaurants## # A tibble: 343 × 6
## x y `n()` x_2 y_2 of_max
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -79.9 40.2 2 -79.9 40.2 0.00844
## 2 -80.0 40.2 1 -79.9 40.2 0.00422
## 3 -79.9 40.2 1 -79.8 40.2 0.00422
## 4 -79.8 40.2 2 -79.8 40.2 0.00844
## 5 -79.9 40.2 1 -79.9 40.3 0.00422
## 6 -79.9 40.2 1 -79.9 40.3 0.00422
## 7 -79.9 40.3 4 -79.9 40.3 0.0169
## 8 -79.9 40.3 30 -79.9 40.3 0.127
## 9 -79.9 40.3 10 -79.8 40.3 0.0422
## 10 -79.8 40.3 1 -79.8 40.3 0.00422
## # … with 333 more rowsOverlay the squares over a map of Pittsburgh using leaflet
library(leaflet)
leaflet() %>%
addTiles() %>%
addRectangles(
restaurants$x,
restaurants$y,
restaurants$x_2,
restaurants$y_2
)Make it fancy… Highlight squares based on number of comparative number of restaurants, and add a popup
fancy <- leaflet() %>%
addTiles() %>%
addRectangles(
restaurants$x,
restaurants$y,
restaurants$x_2,
restaurants$y_2,
fillOpacity = restaurants$of_max,
fillColor = "red",
stroke = FALSE,
popup = paste0("<p>Restaurants: ", restaurants$`n()`,"</p>")
)
fancyKmeans with modeldb
Calculate k-means inside the database using simple_kmeans_db()
service_kmeans <- service %>%
filter(description == restaurant) %>%
simple_kmeans_db(x, y)
head(service_kmeans, 10)## # Source: lazy query [?? x 26]
## # Database: Microsoft SQL Server 12.00.2195[MJB357@pitt.edu@paccmdb/paccmdb]
## k_center k_x k_y id facility_name num street city state zip
## <chr> <dbl> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 center_1 -80.1 40.5 1.07e 4 Bridgeville Da… 413 Bower … Brid… PA 15017
## 2 center_3 -80.0 40.4 1.06e 4 Anthonys Pizza 1820 Browns… Pitt… PA 15210
## 3 center_3 -80.0 40.4 2.00e11 The Pickle Bar… 1301 E Cars… Pitt… PA 15203
## 4 center_2 -79.8 40.4 2.00e11 Cardamone's Pi… 546 Braddo… Brad… PA 15104
## 5 center_3 -80.0 40.4 2.00e11 Pittsburgh New… 3233 Bright… Pitt… PA 15212
## 6 center_1 -80.1 40.5 3.69e 4 Coraopolis Tow… 951 1st Ave Cora… PA 15108
## 7 center_2 -79.8 40.4 5.00e 4 Donut Shack & … 7315 Saltsb… Pitt… PA 15235
## 8 center_2 -79.8 40.4 3.39e 4 Aabco Swansons… 9915 Franks… Pitt… PA 15235
## 9 center_2 -79.8 40.4 3.68e 4 Dumont's Donuts 534 Monong… Glas… PA 15045
## 10 center_3 -80.0 40.4 3.38e 4 Frank & Shirle… 2209 Saw Mi… Pitt… PA 15210
## # … with 16 more variables: municipal <chr>, category_cd <dbl>,
## # description <chr>, p_code <dbl>, fdo <date>, bus_st_date <date>,
## # bus_cl_date <date>, seat_count <dbl>, noroom <dbl>, sq_feet <dbl>,
## # status <chr>, placard_st <dbl>, x <dbl>, y <dbl>, address <chr>, year <dbl>Preview the final SQL statement that the routine settled on
show_query(service_kmeans)## <SQL>
## SELECT "center" AS "k_center", "k_x", "k_y", "id", "facility_name", "num", "street", "city", "state", "zip", "municipal", "category_cd", "description", "p_code", "fdo", "bus_st_date", "bus_cl_date", "seat_count", "noroom", "sq_feet", "status", "placard_st", "x", "y", "address", "year"
## FROM (SELECT "center" AS "k_center", "x" AS "k_x", "y" AS "k_y"
## FROM (SELECT "center", AVG("x") AS "x", AVG("y") AS "y"
## FROM (SELECT "x", "y", "center"
## FROM (SELECT "id", "facility_name", "num", "street", "city", "state", "zip", "municipal", "category_cd", "description", "p_code", "fdo", "bus_st_date", "bus_cl_date", "seat_count", "noroom", "sq_feet", "status", "placard_st", "x", "y", "address", "year", "center_1", "center_2", "center_3", CASE
## WHEN ("center_1" >= "center_1" AND "center_1" < "center_2" AND "center_1" < "center_3") THEN ('center_1')
## WHEN ("center_2" < "center_1" AND "center_2" >= "center_2" AND "center_2" < "center_3") THEN ('center_2')
## WHEN ("center_3" < "center_1" AND "center_3" < "center_2" AND "center_3" >= "center_3") THEN ('center_3')
## END AS "center"
## FROM (SELECT "id", "facility_name", "num", "street", "city", "state", "zip", "municipal", "category_cd", "description", "p_code", "fdo", "bus_st_date", "bus_cl_date", "seat_count", "noroom", "sq_feet", "status", "placard_st", "x", "y", "address", "year", SQRT(((-80.1175574169212 - "x") * (-80.1175574169212 - "x")) + ((40.4811571474466 - "y") * (40.4811571474466 - "y"))) AS "center_1", SQRT(((-79.840383094351 - "x") * (-79.840383094351 - "x")) + ((40.4231194070147 - "y") * (40.4231194070147 - "y"))) AS "center_2", SQRT(((-79.9829052934728 - "x") * (-79.9829052934728 - "x")) + ((40.4354433796054 - "y") * (40.4354433796054 - "y"))) AS "center_3"
## FROM "test_agh_food"
## WHERE ("description" = 'Restaurant without Liquor')) "q01") "q02"
## WHERE (NOT((("center") IS NULL)))) "q03"
## GROUP BY "center") "q04") "LHS"
## RIGHT JOIN (SELECT "id", "facility_name", "num", "street", "city", "state", "zip", "municipal", "category_cd", "description", "p_code", "fdo", "bus_st_date", "bus_cl_date", "seat_count", "noroom", "sq_feet", "status", "placard_st", "x", "y", "address", "year", "center"
## FROM (SELECT "id", "facility_name", "num", "street", "city", "state", "zip", "municipal", "category_cd", "description", "p_code", "fdo", "bus_st_date", "bus_cl_date", "seat_count", "noroom", "sq_feet", "status", "placard_st", "x", "y", "address", "year", "center_1", "center_2", "center_3", CASE
## WHEN ("center_1" >= "center_1" AND "center_1" < "center_2" AND "center_1" < "center_3") THEN ('center_1')
## WHEN ("center_2" < "center_1" AND "center_2" >= "center_2" AND "center_2" < "center_3") THEN ('center_2')
## WHEN ("center_3" < "center_1" AND "center_3" < "center_2" AND "center_3" >= "center_3") THEN ('center_3')
## END AS "center"
## FROM (SELECT "id", "facility_name", "num", "street", "city", "state", "zip", "municipal", "category_cd", "description", "p_code", "fdo", "bus_st_date", "bus_cl_date", "seat_count", "noroom", "sq_feet", "status", "placard_st", "x", "y", "address", "year", SQRT(((-80.1175574169212 - "x") * (-80.1175574169212 - "x")) + ((40.4811571474466 - "y") * (40.4811571474466 - "y"))) AS "center_1", SQRT(((-79.840383094351 - "x") * (-79.840383094351 - "x")) + ((40.4231194070147 - "y") * (40.4231194070147 - "y"))) AS "center_2", SQRT(((-79.9829052934728 - "x") * (-79.9829052934728 - "x")) + ((40.4354433796054 - "y") * (40.4354433796054 - "y"))) AS "center_3"
## FROM "test_agh_food"
## WHERE ("description" = 'Restaurant without Liquor')) "q01") "q02"
## WHERE (NOT((("center") IS NULL)))) "RHS"
## ON ("LHS"."k_center" = "RHS"."center")Create a local table with the corners of each center
cr <- service_kmeans %>%
group_by(k_x, k_y) %>%
summarise(
restaurants = n(),
lat1 = min(y, na.rm = TRUE),
lat2 = max(y, na.rm = TRUE),
lon1 = min(x, na.rm = TRUE),
lon2 = max(x, na.rm = TRUE),
) %>%
collect()## `summarise()` has grouped output by 'k_x'. You can override using the `.groups` argument.cr## # A tibble: 3 × 7
## # Groups: k_x [3]
## k_x k_y restaurants lat1 lat2 lon1 lon2
## <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
## 1 -79.8 40.4 1376 40.2 40.7 -79.9 -79.7
## 2 -80.1 40.5 849 40.3 40.7 -80.3 -80.0
## 3 -80.0 40.4 2753 40.2 40.7 -80.1 -79.9Plot against a map
leaflet() %>%
addTiles() %>%
addRectangles(cr$lon1, cr$lat1, cr$lon2, cr$lat2) Overlay the raster boxes
leaflet() %>%
addTiles() %>%
addRectangles(cr$lon1, cr$lat1, cr$lon2, cr$lat2) %>%
addRectangles(
restaurants$x,
restaurants$y,
restaurants$x_2,
restaurants$y_2,
fillColor = "red",
stroke = FALSE,
popup = paste0("<p>Restaurants: ", restaurants$`n()`,"</p>")
) Bringing it all together
Use purrr to run db_compute_raster() for each center in the kmeans results
library(purrr)
cents <- c("center_1", "center_2", "center_3")
by_centers <- cents %>%
map_df(~{
bc <- service_kmeans %>%
filter(k_center == .x) %>%
db_compute_raster2(x, y, resolution = 20) %>%
mutate(
of_max = `n()` / max(`n()`),
center = .x
)
})Change the color of each center
bc <- by_centers %>%
mutate(color = case_when(
center == "center_1" ~ "red",
center == "center_2" ~ "blue",
center == "center_3" ~ "green",
TRUE ~ "white"
))
bc## # A tibble: 424 × 8
## x y `n()` x_2 y_2 of_max center color
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 -80.1 40.3 1 -80.1 40.3 0.0196 center_1 red
## 2 -80.1 40.3 1 -80.1 40.3 0.0196 center_1 red
## 3 -80.1 40.3 5 -80.1 40.4 0.0980 center_1 red
## 4 -80.1 40.3 2 -80.1 40.4 0.0392 center_1 red
## 5 -80.2 40.4 10 -80.1 40.4 0.196 center_1 red
## 6 -80.1 40.4 35 -80.1 40.4 0.686 center_1 red
## 7 -80.1 40.4 8 -80.1 40.4 0.157 center_1 red
## 8 -80.1 40.4 9 -80.1 40.4 0.176 center_1 red
## 9 -80.2 40.4 2 -80.2 40.4 0.0392 center_1 red
## 10 -80.1 40.4 10 -80.1 40.4 0.196 center_1 red
## # … with 414 more rowsMap the two layers together
leaflet() %>%
addTiles() %>%
addRectangles(cr$lon1, cr$lat1, cr$lon2, cr$lat2) %>%
addRectangles(
bc$x,
bc$y,
bc$x_2,
bc$y_2,
fillColor = bc$color,
fillOpacity = 0.5,
stroke = FALSE,
popup = paste0("<p>Restaurants: ", bc$`n()`,"</p>")
)dbDisconnect(con)