In this file we will continue practicing how to query online databases with SQL. We will use our BigQuery setup from the previous document with two public databases to demonstrate queries for large datasets. We will also practice handling spatial data, and making maps.

Loading packages:

```{r}
# For bigquery and general munging
library("DBI")
library("bigrquery")
library("tidyverse")
# And some packages for spatial data and mapping (see also "ggmaps", "sf", and "rgdal")
library("sf")
library("tmap")
library("classInt")
```

## London LSOA Crime statistics

We are first going to explore a dataset of London crime statistics. The data includes the number of crimes at two different geographic levels of London (LSOA and borough) by year, according to crime type. First, we have to change our project name such that the `bigrquery` package knows that it has to request public datasets (previously we were requesting our own datasets).

```{r}
# Project name (now accessing the public datasets)
project_name <- "bigquery-public-data"

# Billing (main project ID)
billing_info <- "tba"

# Dataset name (London crime database within the public databases)
dataset_name <- "london_crime"
```

1. Creating the database connection object, and exploring the dataset by field.

```{r}
db <- dbConnect(
  bigrquery::bigquery(),
  project = project_name,
  dataset = dataset_name,
  billing = billing_info
)

# See what tables are available in this dataset
dbListTables(db)

# See the fields in the table of interest
dbListFields(db, "crime_by_lsoa")
```

2. Let's count how many rows the table contains. Remember, if the code yields a Google credentials error, once copy-paste it into the R console directly which should then open the browser window for authentication -- alternatively try and reinstall `rlang`.

```{r}


dbGetQuery(db, "SELECT COUNT(*) FROM crime_by_lsoa")
```

3. We can also extract crimes per year

```{r}
dbGetQuery(db, "SELECT year, SUM(value) AS count_crime FROM crime_by_lsoa
           GROUP BY year
           ORDER BY year")
```

4. And the same but by year **and** borough

```{r}
data_plot_1 <- dbGetQuery(db, "SELECT year, borough, SUM(value) AS count_crime
            FROM crime_by_lsoa
            GROUP BY year, borough")
View(data_plot_1)
```

```{r fig.width=7, fig.height=4, echo=FALSE}
data_plot_1 |> 
  group_by(borough) |>
  mutate(ratio = count_crime/count_crime[year == 2008]) |>
  ggplot(aes(x = year, y = ratio, color = borough)) + 
    geom_line() + #add lines
    labs(title = "Crime in London's boroughs, over time", x = "Year", y = "Crime in Year = T / Crime in 2008") + #label the plot
    theme_minimal() +
    theme(legend.key.size = unit(0.01, "npc")) + #relatively rescale the size of the legend 
    labs(color = "Borough") #re-label the legend
```

5. Crime by year and category

```{r}
data_plot_2 <- dbGetQuery(db, "SELECT year, major_category,
            SUM(value) AS count_crime FROM crime_by_lsoa
            GROUP BY year, major_category")
#View(data_plot_2)
```

```{r}
ggplot(data_plot_2) + aes(x = year, y = count_crime, colour = major_category) +
  geom_line() +
  theme_minimal() +
  labs(title = "Crime in London by major category, over time", y = "Crime count", x = "Year") + #label the plot
  labs(color = "Major Category") #re-label the legend
```

6. Finally, let's make a map that depicts crime in 2016 across the city. First, let's call in the crime data for 2016 at the LSOA level, which is nice and granular so will be useful for a map.

```{r}
# Call crime at the Lower Layer Super Output Area level
crime_per_lsoa_code <- dbGetQuery(db, "SELECT lsoa_code, value AS count_crime
            FROM crime_by_lsoa
            WHERE year = 2016
            GROUP BY lsoa_code")
crime_per_lsoa_code
```

Now we need some way to **map** this data. To do that, we need an object that has what is called a **geometry** -- a **mapping** of our data onto a two-dimensional representation over space. These files are typically called `shapefiles`, and are typically stored as `.shp` files that are connected to a series of related files of the same name but different extension (e.g. `.dbf` files, which are dBASE tables that store attributes/features/variables). We will download a shapefile for London's LSOAs from the London Government website, and then extract that `.zip` file (which includes a folder structure and many other files) into our data folder, creating a new sub-folder called `shapefiles`.  

```{r}
# Download a shapefile from the London Government website:
url_dl <- "https://data.london.gov.uk/download/statistical-gis-boundary-files-london/9ba8c833-6370-4b11-abdc-314aa020d5e0/statistical-gis-boundaries-london.zip"
file_name_dl <- "statistical-gis-boundaries-london.zip"
file_path_dl <- "data/"

download.file(url_dl, paste0(file_path_dl, file_name_dl), mode = "wb")

# Unzip file: 
unzip(paste0(file_path_dl, file_name_dl),exdir=paste0(file_path_dl,"shapefiles"))

# Read the shapefile:
shp <- st_read("data/shapefiles/statistical-gis-boundaries-london/ESRI/LSOA_2011_London_gen_MHW.shp")

# Let's explore the shapefile a little: 
class(shp)
names(shp)
head(shp)
```

Now that we have the shapefile downloaded, let's combine the shapefile and the data we queried earlier, by the LSOA code. 

```{r}
# Merging our crime data into the shapefile:
shp_to_plot <- shp |>
  left_join(crime_per_lsoa_code, by = join_by(LSOA11CD == lsoa_code))

dim(shp_to_plot)
head(shp_to_plot)
```

Now we can start to map things:

```{r}
# First, let's just look at the geometry
plot(st_geometry(shp_to_plot))

# Now, a *very* simple plot of crime count... does this look good?
plot(shp_to_plot["count_crime"])

# We can also make a simple plot with ggplot2 (see 'ggmap' for more advanced usage, though you will likely then need to use 'sp' instead of 'sf')
map <- ggplot(shp_to_plot) +
  geom_sf(aes(fill=count_crime)) +
  theme_minimal()

map
```

This was fine, but we can approach this in a slightly more bespoke way, using [Tmap](https://r-tmap.github.io/tmap-book/).

```{r}
# Alternatively, we can use the 'tmaps' package, which is a powerful package built for thematic mapping 
map <- tm_shape(shp_to_plot) +
  tm_polygons("count_crime")

map

# Let's improve the map by creating more visual variation with our bins: 
breaks <- classInt::classIntervals(shp_to_plot$count_crime, n = 5, style = "jenks", dig.lab=10)
breaks

shp_to_plot <- shp_to_plot |>
  mutate(count_crime_cut = cut(count_crime, breaks$brks, include.lowest = TRUE, dig.lab=10))

map <- tm_shape(shp_to_plot) +
  tm_polygons("count_crime_cut", palette = "Reds", title = "Number of Crimes") +
  tm_layout(main.title = "Crime in London, 2016", legend.outside = TRUE, frame = FALSE)
  
map
tmap_save(map, "london_crime.png", dpi = 300)

# Finally, we can create a dynamic "live" map using tmaps:
tmap_options(check.and.fix = TRUE) 
tmap_mode("view")
map

tmap_save(map, "london_crime.html")
tmap_mode("plot")
```

## NYC Bicycle Hire

The second dataset contains a table about Citi Bike trips (NYC's bike-sharing service) since Citi Bike launched in September 2013.

First, we need to update the dataset name and connect again. Then we should explore the dataset, and see what tables exist.

```{r}
dataset_name <- "new_york"

db <- dbConnect(
  bigrquery::bigquery(),
  project = project_name,
  dataset = dataset_name,
  billing = billing_info
)

# Check what tables exist
dbListTables(db)

# Explore the field in the table of interest
dbListFields(db, "citibike_trips")
```

1. First, count how many rows the table of interest contains:

```{r}
# Number of trips in the table of interest
dbGetQuery(db, "SELECT COUNT(*) FROM citibike_trips")
```

```{r}
# First 10 entries
dbGetQuery(db, "SELECT * FROM citibike_trips LIMIT 10")
```

2. Which are the 10 most popular stations in terms of how many trips started there?

```{r}
dbGetQuery(db, "SELECT start_station_name, COUNT(*) AS count_start
  FROM citibike_trips
  GROUP BY start_station_name
  ORDER BY count_start DESC
  LIMIT 10")
```

3. What is the average trip duration of a CitiBike trip in NYC? For advanced keywords like these, see the detailed [documentation](https://cloud.google.com/bigquery/docs/reference/standard-sql/timestamp_functions).

```{r}
dbGetQuery(db, "SELECT AVG( TIMESTAMP_DIFF(stoptime, starttime, MINUTE) )
           AS duration_minutes
           FROM citibike_trips")
```

4. What is the average trip duration based on the hour of the day when they start?

```{r}
dbGetQuery(db, "SELECT EXTRACT(HOUR FROM starttime) as hour_of_day,
  AVG( TIMESTAMP_DIFF(stoptime, starttime, MINUTE) ) AS duration_minutes
  FROM citibike_trips
  GROUP BY hour_of_day
  ORDER BY hour_of_day")
```

5. What is the average trip duration grouped by the gender column?

```{r}
dbGetQuery(db, "SELECT gender, AVG( TIMESTAMP_DIFF(stoptime, starttime, MINUTE) )
  AS duration_minutes, COUNT(*) as trip_count,
  FROM citibike_trips
  GROUP BY gender")
```

6. What is the average distance of a trip?

```{r}
dbGetQuery(db, "SELECT AVG( (ABS(start_station_latitude-end_station_latitude) +
      ABS(start_station_longitude-end_station_longitude)) * 111) AS avg_distance_km
  FROM citibike_trips")
```

7. What is the average distance of a trip grouped by the gender column?

```{r}
dbGetQuery(db, "SELECT gender, AVG( (ABS(start_station_latitude-end_station_latitude) +
      ABS(start_station_longitude-end_station_longitude)) * 111) AS avg_distance_km
  FROM citibike_trips
  GROUP BY gender")
```





















