Cran Contributed Packages

library(maptools)
library(rgdal)
library(sp)
library(sf)
library(tigris)
library(hrbrthemes)
library(tidyverse)
require(spatialEco)
library(maptools)
###library(GISTools)
library(rgdal)


## March Crash
setwd("D:/From Syncplicity/Bulk_Papers/03142019/covid/NY_TrafficCrash/covid-19-data-master")
march <- read.csv("NY_March_CrashesLL03.csv")
dim(march)
## [1] 117633      5
names(march)
## [1] "ID"             "Year"           "LATITUDE"       "LONGITUDE"     
## [5] "ON.STREET.NAME"
coordinates(march)=~LONGITUDE+LATITUDE
proj4string(march)<- CRS("+proj=longlat +datum=WGS84")
march1 <-spTransform(march,CRS("+proj=longlat"))


### NY TaxiZone
st_read("D:/From Syncplicity/Bulk_Papers/03142019/covid/NY_TrafficCrash/zones.geojson") %>%
  st_set_crs(4326) %>%
  st_transform("+proj=longlat +datum=WGS84 +no_defs") -> ny
## Reading layer `zones' from data source `D:\From Syncplicity\Bulk_Papers\03142019\covid\NY_TrafficCrash\zones.geojson' using driver `GeoJSON'
## Simple feature collection with 263 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
## CRS:            4326
plot(ny)

pts_sf = st_as_sf(march1)
p_sf = st_as_sf(ny)

st_crs(p_sf)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]
p_sf1= st_transform(p_sf, st_crs(pts_sf))

##https://mgimond.github.io/Spatial/index.html
## https://cengel.github.io/R-spatial/spatialops.html
nn= p_sf1 %>% 
  st_join(pts_sf) %>%
  group_by(zone, Year) %>% 
  summarize(count = n()) 

nn1=subset(nn, Year==2020)
nn2=subset(nn, Year==2019)
colnames(nn1)[3] <- "2020 Crashes"
p1= plot(nn1["2020 Crashes"])

colnames(nn2)[3] <- "2019 Crashes"
p2= plot(nn2["2019 Crashes"])

#par(mfrow=c(1,2))
#plot(nn1["2020 Crashes"])
#plot(nn2["2019 Crashes"])
library(tmap)

tm_shape(nn) +
  tm_polygons("count", palette = "RdYlBu") +
  tm_facets(by = "Year")

nn1= subset(nn, Year> 2016 & Year < 2021)

tm_shape(nn1) +
  tm_polygons("count", palette = "RdYlBu") +
  tm_facets(by = "Year")