Shapefile source: https://vangdata.carto.com/tables/shapefiles_barcelona_distrito/public

#install sf from github using:
#install.packages("devtools")
#devtools::install_github("r-spatial/sf")
#everything else can be installed normally
library(ggplot2)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggmap)
## Warning: package 'ggmap' was built under R version 4.2.3
## ℹ Google's Terms of Service: ]8;;https://mapsplatform.google.com<https://mapsplatform.google.com>]8;;
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
#load in the base map - set this to your own file path
setwd("~/Binghamton/geog380/barcelona_map")
barcelona <- st_read("shapefiles_barcelona_distrito.shp")
## Reading layer `shapefiles_barcelona_distrito' from data source 
##   `C:\Users\mhaller\Documents\Binghamton\geog380\barcelona_map\shapefiles_barcelona_distrito.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2.052451 ymin: 41.31692 xmax: 2.229203 ymax: 41.46825
## Geodetic CRS:  WGS 84
#select a projection
#don't worry too much about the code - this is just a projection designed for Spain
barcelona <- barcelona %>% st_transform("EPSG:2062")

#plot using ggplot
ggplot()+
  geom_sf(barcelona, mapping=aes(), fill = "white")

#load in your points
setwd("~/Binghamton/geog380")
accidents <- read.csv("accidents_2017.csv")

#convert to sf object
accidents1 <- accidents %>% 
  #I need to transform my lat/long coordinates into an sf object 
  #4326 is the crs used in gps - it is defined in lat/long coords
  #we need to first convert to an sf object using this crs
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>% 
  #now we can correctly convert to the same crs as our barcelona map
  st_transform(crs = st_crs(barcelona))

#let's plot both layers
ggplot()+
  geom_sf(barcelona, mapping=aes(), fill = "gray80")+
  geom_sf(accidents1, mapping=aes(color = District.Name))+
  theme_minimal()