R Notebook

install.packages(c(tidyverse, sf, tmap, spatstat, maptools, ggplot2, sqldf)) install.packages(“spatstat.geom”, dependencies = TRUE) install.packages(“leaflet.providers”) install.packages(“markdown”, dependencies = TRUE) install.packages(“polyclip”, dependencies = TRUE) install.packages(“ggmap”, dependencies = TRUE)

library(tidyverse) 
library(sf)
library(tmap)
library(spatstat)

Loading required package: spatstat.data
Loading required package: spatstat.geom
Registered S3 method overwritten by 'spatstat.geom':
  method     from
  print.boxx cli 
spatstat.geom 2.2-2
Loading required package: spatstat.core
Loading required package: nlme

Attaching package: 㤼㸱nlme㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    collapse

Loading required package: rpart
spatstat.core 2.3-0
Loading required package: spatstat.linnet
spatstat.linnet 2.3-0

spatstat 2.2-0       (nickname: 㤼㸱That's not important right now㤼㸲) 
For an introduction to spatstat, type 㤼㸱beginner㤼㸲 

library(maptools)

Loading required package: sp
Checking rgeos availability: FALSE
    Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
    which has a restricted licence. It is disabled by default;
    to enable gpclib, type gpclibPermit()

library(sqldf)

Loading required package: gsubfn
Loading required package: proto
Loading required package: RSQLite

library(ggplot2)
library(spatstat.geom)
library(leaflet.providers)
library(markdown)
library(polyclip)
library(ggmap)

Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
Please cite ggmap if you use it! See citation("ggmap") for details.

######## PART 3 SEGMENTATION AND MAPPING ######### 

# Sample for just the morning 
morning <- rdfClean %>% filter(hour %in% c('07', '08', '09'))

# Find most important stations
morningStations <- morning  %>%
  group_by(`start.station.id`) %>%
  summarise(lat=as.numeric(`start.station.latitude`[1]),
            lon=as.numeric(`start.station.longitude`[1]),
            name=`start.station.name`[1],
            trips=n())

## Take a look
morningStations %>%
  arrange(desc(trips))

tmap_mode('view')

tmap mode set to interactive viewing

tm_shape(morningStationsGeo %>% 
           filter(!is.na(trips)) %>% 
           mutate(trips = as.numeric(trips))) + tm_dots(col = "trips", style="quantile", id = "name")

## Read SF for a file of boroughs
boroughs <- st_read("borough.shp")

Reading layer `borough' from data source 
  `C:\Users\dpoddar3\Downloads\Assignment__Citibikes_export\ForClass\borough.shp' using driver `ESRI Shapefile'
Simple feature collection with 5 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
Geodetic CRS:  WGS84(DD)

boroughs_prj <- st_transform(boroughs, crs = 2263)

# Density map (This method doesn't allow us to see the background map)
df.size.proj <- st_transform(morningStationsGeo, crs = 2263)
p.sp <- as(df.size.proj, "Spatial")
p.ppp <- as(p.sp, "ppp")
kernelDensityMorning <- (density(p.ppp, sigma = 1000))

##Plot the two together
plot(kernelDensityMorning)
plot(boroughs_prj$geometry, add=TRUE)

##CHALLENGES FOR YOUR DENSITY MAP
##Extra: Clip the raster and clip the vector as well. (only plot the intersection of the two?) 
##Extra: Make the lines white (hint: it is now a polygon!)
##Extra: Add the top 3 (or 5) points and label them.

# Heatmap in ggplot
library(ggplot2)
library(ggmap)

df.size.degree <- st_transform(morningStationsGeo, crs = 4326) %>% 
  mutate(lat = st_coordinates(.)[,2],
         lon = st_coordinates(.)[,1])

##Get borders
height <- max(st_coordinates(df.size.degree)[,2]) - min(st_coordinates(df.size.degree)[,2])
width <- max(st_coordinates(df.size.degree)[,1]) - min(st_coordinates(df.size.degree)[,1])
  
sac_borders <- c(bottom = min(st_coordinates(df.size.degree)[,2]) - 0.1*height,
                top = max(st_coordinates(df.size.degree)[,2]) + 0.1*height,
                left = min(st_coordinates(df.size.degree)[,1]) - 0.1*width,
                right = max(st_coordinates(df.size.degree)[,1]) + 0.1*width)

map <- get_stamenmap(sac_borders, zoom = 12, maptype = "toner-lite")

Source : http://tile.stamen.com/toner-lite/12/1205/1538.png
Source : http://tile.stamen.com/toner-lite/12/1206/1538.png
Source : http://tile.stamen.com/toner-lite/12/1205/1539.png
Source : http://tile.stamen.com/toner-lite/12/1206/1539.png
Source : http://tile.stamen.com/toner-lite/12/1205/1540.png
Source : http://tile.stamen.com/toner-lite/12/1206/1540.png
Source : http://tile.stamen.com/toner-lite/12/1205/1541.png
Source : http://tile.stamen.com/toner-lite/12/1206/1541.png

ggmap(map) +
  geom_density2d(data = df.size.degree, aes(x = lon, y = lat), size = 0.3) +
  stat_density2d(data = df.size.degree, aes(x = lon, y = lat, fill = ..level.., alpha = ..level..), 
                 geom = "polygon", bins = 10) + 
  scale_fill_gradient(low = "green", high = "red") +
  scale_alpha(range = c(0, 0.4), guide = FALSE) #Plot +

 coord_cartesian(xlim = c(-74.02, -95.7), 
                  ylim = c(29.5, 30.1))

<ggproto object: Class CoordCartesian, Coord, gg>
    aspect: function
    backtransform_range: function
    clip: on
    default: FALSE
    distance: function
    expand: TRUE
    is_free: function
    is_linear: function
    labels: function
    limits: list
    modify_scales: function
    range: function
    render_axis_h: function
    render_axis_v: function
    render_bg: function
    render_fg: function
    setup_data: function
    setup_layout: function
    setup_panel_guides: function
    setup_panel_params: function
    setup_params: function
    train_panel_guides: function
    transform: function
    super:  <ggproto object: Class CoordCartesian, Coord, gg>