Demonstrating Advanced Map Knowledge through Philidelphia Crime Data

1 Case Study

For the case study we will be looking at Philidelphia Neighborhood Shootings. We will be getting this data in the form of a geojson file.

phillyNeighborShooting  <- na.omit(st_read("https://pengdsci.github.io/STA553VIZ/w08/PhillyShootings.geojson"))
phillyNeighbor  <- st_read("https://pengdsci.github.io/STA553VIZ/w08/Neighborhoods_Philadelphia.geojson")
philly  <- st_read("https://pengdsci.github.io/STA553VIZ/w08/PhillyNeighborhood-blocks.geojson") # block level data
SubCrime = phillyNeighborShooting[,c("dc_key","fatal","location","year")]
aggregateCrime = aggregate(SubCrime$fatal, by=list(SubCrime$location,SubCrime$year), FUN=length)
## Longitude and latitude of zip code center
ZipLon = aggregate(as.numeric(phillyNeighborShooting$point_x), by=list(phillyNeighborShooting$location), FUN=mean)
ZipLat = aggregate(as.numeric(phillyNeighborShooting$point_y), by=list(phillyNeighborShooting$location), FUN=mean)
ZipLonLat = merge(ZipLon, ZipLat, by = "Group.1")
names(ZipLonLat) = c("location", "lng", "lat")
##  Zip code crime
ZipCrimeTable = table(phillyNeighborShooting$fatal,  phillyNeighborShooting$location)
ZipCrime = data.frame(location=names(table(phillyNeighborShooting$location)), 
                      fatal = table(phillyNeighborShooting$fatal, phillyNeighborShooting$location)[1,],
                      nonfatal = table(phillyNeighborShooting$fatal, phillyNeighborShooting$location)[2,],
                      total.crime = table(phillyNeighborShooting$location) 
                      )
ZipCrime = ZipCrime[, c("location", "fatal", "nonfatal", "total.crime.Freq")]
colnames(ZipCrime) = c("location", "fatal", "nonfatal", "total.crime")
ZipCrime = merge(ZipLonLat, ZipCrime, by = "location")

2 Creating the plots to be turned into images

2.1 Wound Histogram

plot_ly(
  data = phillyNeighborShooting,
   x = ~ tolower(wound),
   type = "histogram",
   nbinsx = 10, 
   name = "Area Wounded",
   alpha = .5,
   marker = list(line = list(color = "darkgray",  width = 2)) ) %>%
   ## Adding additional histograms and stacking them
    add_histogram(x = ~tolower(wound),
                name = "Area Wounded", nbinsx = 10, alpha = 0.5,
                marker = list(line = list(color = "darkgray",  width = 2))) %>%
    layout(barmode = "overlay",
         title = "Histogram of Wounded Areas",
         xaxis = list(title = "Area Wounded",
                      zeroline = TRUE),
         yaxis = list(title = "Count",
                      zeroline =TRUE),
               ## margin of the plot
      margin = list(
              b = 50,
              l = 100,
              t = 120,
              r = 50
      ))

2.2 Box Plot of Fatalities for each Race

barplotdata = aggregate(phillyNeighborShooting[,21], by = list(phillyNeighborShooting$race), FUN = sum)

plot_ly(
  data = barplotdata,
   x = ~Group.1,
   y = ~fatal,
   type = "bar",
   name = "fatal",
   ## graphic size
   width = 700,
   height = 400) %>%
    layout( yaxis = list(title ="Number of Fatalities"),
            xaxis = list(title = "Race"),
            title = "Fatalaty Shootings by Race"
       )

2.3 Density of Sex

# define three densities
male.which <- phillyNeighborShooting[which(phillyNeighborShooting$sex == "M"),]
male <- density(male.which$fatal)
female.which <- phillyNeighborShooting[which(phillyNeighborShooting$sex == "F"),]
female <- density(female.which$fatal)
# plot density curves
fig <- plot_ly(x = ~male$x, 
               y = ~male$y, 
               type = 'scatter', #A character string specifying the trace type
               mode = 'lines', 
               name = 'males', 
               fill = 'tozeroy')  %>% 
           # adding more density curves
       add_trace(x = ~female$x, 
                 y = ~female$y, 
                 name = 'female', 
                 fill = 'tozeroy')  %>%   
       layout(xaxis = list(title = 'Sex'),
              yaxis = list(title = 'Density'))
fig

3 Final Map

## Images to be plot on the map
fatalitybyrace = "https://dannyoc99.github.io/DOCAL/img/FatalitybyRace.png"
SexFatality = "https://dannyoc99.github.io/DOCAL/img/SexFatality.png"
WoundHistogram = "https://dannyoc99.github.io/DOCAL/img/WoundHistogram.png"
fata = st_as_sf(data.frame(x = -75.3877, y = 39.9168),
                coords = c("x", "y"),
                crs = 4326)
SexFatalityloc = st_as_sf(data.frame(x = -75.3677, y = 39.9168),
                coords = c("x", "y"),
                crs = 4326)
WoundHistogramloc = st_as_sf(data.frame(x = -75.3477, y = 39.9168),
                coords = c("x", "y"),
                crs = 4326)


##################################################
## Crime distribution by zip code - bubble plot
pal0 <- colorNumeric(palette = viridis(256, option = "B"), domain = range(ZipCrime$total.crime))
zipfig <- leaflet() %>%
  setView(lng=-75.1527, lat=39.9707, zoom = 11) %>%
  addProviderTiles(providers$Esri.WorldGrayCanvas) %>% 
  ## plot the neighborhood boundary of Philly
  addPolygons(data = phillyNeighbor,
              color = 'skyblue',
              weight = 1)  %>%
  ##
  addCircleMarkers(data = ZipCrime,
                   radius = ~((total.crime)^(1/3)),
                   stroke = FALSE, 
                   fillOpacity = 0.5,
                   popup = ~paste('Location:', location, 
                      '<br>Total Crime:', total.crime, 
                      '<br>Fatal Crime:', fatal,
                      '<br>Nonfatal Crime:', nonfatal)) 
##################################################
## popup interactive graphs with plotly
fl = tempfile(fileext = ".html")
saveWidget(zipfig, file = fl)
###
leaflet() %>%
  setView(lng=-75.15092, lat=40.00995, zoom = 11) %>%
  #addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  ##
  addProviderTiles(providers$CartoDB.DarkMatter, group="Dark") %>%
  addProviderTiles(providers$CartoDB.DarkMatterNoLabels, group="DarkLabel") %>%  
  addProviderTiles(providers$Esri.NatGeoWorldMap, group="Esri") %>%
  #addTiles(providers$CartoDB.PositronNoLabels) %>%
  ## mini reference map
  addMiniMap() %>%
  ## neighborhood boundary
  addPolygons(data = phillyNeighbor,
              color = 'skyblue',
              weight = 1)  %>%
  ## plot information on the map
  addCircleMarkers(data = phillyNeighborShooting,
                   radius = ~ifelse(fatal == "1", 5, 3),
                   stroke = FALSE, 
                   fillOpacity = 0.5,
                   popup = ~popupTable(phillyNeighborShooting),
                   clusterOptions = markerClusterOptions(maxClusterRadius = 40)) %>%

  # Trend of crimes over the years
  addCircleMarkers(data = fata, 
                   color = "red",
                   weight = 2,
                   label = "fatality by race",
                   stroke = FALSE, 
                   fillOpacity = 0.95,
                   group = "fata") %>%
  addPopupImages(fatalitybyrace, 
                  width = 500,
                  height = 350,
                  tooltip = FALSE,
                  group = "fata") %>%
 # age distribution within the two types of crimes
 addCircleMarkers(data = SexFatalityloc, 
                   color = "black",
                   weight = 2,
                   label = "Density of Fatalities by Sex",
                   stroke = FALSE, 
                   fillOpacity = 0.95,
                   group = "SexFatalityloc") %>%
  addPopupImages(SexFatality, 
                  width = 500,
                  height = 320,
                  tooltip = FALSE,
                  group = "SexFatalityloc") %>%
  addCircleMarkers(data = WoundHistogramloc, 
                   color = "blue",
                   weight = 2,
                   label = "Wound Location Histogram",
                   stroke = FALSE, 
                   fillOpacity = 0.95,
                   group = "WoundHistogramloc") %>%
  addPopupImages(WoundHistogram, 
                  width = 500,
                  height = 400,
                  tooltip = FALSE,
                  group = "WoundHistogramloc") %>% 
  ##
  browsable()

As a case study I decided to graph the shootings by location on a map of the Philidelphia Area. By looking at the map we can see that a crimes are generally in the Center-North Philidelphia Area. Given our data, I was curious to see whether or not the information about these shootings could be expanded on so I made 3 graphs to get more information out of our map besides where these crimes were taking place. Our first graph, fatality by race shows that a large majority of our fatalities are African American shootings. The next graph, density of fatal shootings by sex shows that while females have a lower density, their density curves are wider in both the fatal and non-fatal shootings. Our final graph, a histogram, shows the count of where people are getting shot. The histogram shows that getting shot in multiple locations is the most commmon, but if you were to limit the shooting to only one area, the leg is the most common place to get shot.

4 Chloreopath Data Preparation

democrat2020<-read_csv(file ="https://dannyoc99.github.io/DOCAL/dataset/countypresidential_election_2000-2020.csv")%>%filter(year==2020)%>%filter(party=="DEMOCRAT")%>%subset(select=c(state_po, county_name, county_fips, candidatevotes))%>%rename(dem_candidate_votes=candidatevotes)
rep2020<-read_csv(file ="https://dannyoc99.github.io/DOCAL/dataset/countypresidential_election_2000-2020.csv")%>%filter(year==2020)%>%filter(party=="REPUBLICAN")%>%subset(select=c(county_fips, candidatevotes))%>%rename(rep_candidate_votes=candidatevotes)
election2020<-inner_join(democrat2020, rep2020, by = c("county_fips" = "county_fips")) 
fips_to_geo<-read_csv(file ="https://dannyoc99.github.io/DOCAL/dataset/fips2geocode.csv")
all_data<-inner_join(election2020, fips_to_geo, by = c("county_fips" = "fips"))
all_data$rep_to_dem<-with(all_data,rep_candidate_votes-dem_candidate_votes)

To look at the 2020 election we will limit our election data to 2020, only look at democrat vs republican candidates. We will also bring in Geo information via an inner join

5 Chloreopath Map of the 2020 Election

#library(plotly)
#install.packages('rjson')
#library("RColorBrewer")  # brewer.pal.info for list of color scales

url <- 'https://raw.githubusercontent.com/plotly/datasets/master/geojson-counties-fips.json'  
counties <- rjson::fromJSON(file=url)  

df=all_data
g <- list(
      scope = 'usa',
      projection = list(type = 'albers usa'),
      showlakes = TRUE,
      lakecolor = toRGB('white')
    )
###
fig <- plot_ly()  %>% 
  add_trace( type = "choropleth",
          geojson = counties,
        locations = df$fips,
                z = df$rep_to_dem,
       colorscale = "RdBu",  
             zmin = 0,
             zmax = 30,
             text = df$name, # hover mesg
           marker = list(line=list(width=0.2))
          )   %>% 
  colorbar(title = "Party (Rep to Dem)",
           colorscale='Viridis')  %>% 
  layout( ### Title 
          title =list(text = "US 2020 Election Choropleth", 
                          font = list(family = "Times New Roman",  # HTML font family  
                                        size = 18,
                                       color = "red")), 
          geo = g)
## The actual HTML page does shoe colorbar correctly, we hide the legend.
hide_legend(fig)

I could not get the graph to display exactly correctly, but the code behind this looks correct, so I thought it may have something to do with the package and thought it worth submitting.