A repro example to get data, and aggregate points into polygons over a list with purrr::map and then animate with ggplot

sf point_in_poly.r

library(corrplot)
library(viridis)
library(stargazer)
library(tidyverse)
library(dplyr)
library(sf)
library(tigris)
library(ggplot2)
library(rgdal)
library(maptools)
library(lubridate)
library(gganimate)
# library(gapminder)
library(transformr)
library(tweenr)
# library(janeaustenr)
library(gifski)
library(magick)
library(rgdal)
# library(QuantPsyc)
library(caret)
library(spatstat)
library(spdep)
library(FNN)
library(grid)
library(gridExtra)
library(purrr)
# library(rCarto)
library(mapedit)

options(scipen = 999)

####Define universal map theme####
mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 14,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  )
}


#Load in Philadelphia Neighborhood Basemap
neighborhoods <- 
  st_read("https://raw.githubusercontent.com/azavea/geo-data/master/Neighborhoods_Philadelphia/Neighborhoods_Philadelphia.geojson") %>%
  st_transform(crs=3365) %>%
  dplyr::select(District = mapname)

tracts <- 
  st_read("http://data.phl.opendata.arcgis.com/datasets/8bc0786524a4486bb3cf0f9862ad0fbf_0.geojson") %>%
  st_transform(crs=3365) %>%
  dplyr::select(District = GEOID10)

bothBoundaries <- rbind(mutate(neighborhoods, Legend = "Neighborhoods"), 
                        mutate(tracts, Legend = "Tracts"))

ggplot() +
  geom_sf(data = bothBoundaries) +
  facet_wrap(~Legend) +
  labs(title = "Local Boundaries, Philadelphia") +
  mapTheme()

phlBoundary <- 
  st_read("http://data.phl.opendata.arcgis.com/datasets/405ec3da942d4e20869d4e1449a2be48_0.geojson") %>%
  st_transform(crs=3365)

####Make Fishnet####
fishnet <- 
  st_make_grid(phlBoundary, cellsize = 1200) %>%
  st_sf()

fishnet <- 
  fishnet[phlBoundary,] %>%
  mutate(uniqueID = rownames(.)) %>%
  dplyr::select(uniqueID)

ggplot() +
  geom_sf(data=fishnet) +
  labs(title = "Fishnet in Philadelphia") +
  mapTheme()

####Load New Construction Permits####
NewCon <-  st_read("https://phl.carto.com/api/v2/sql?q=SELECT+*+FROM+li_permits+WHERE+typeofwork+=+%27NEWCON%27&filename=li_permits&format=geojson&skipfields=cartodb_id") %>%  
  filter(status != "REVOKED") %>% 
  filter(ownername != "PHILADELPHIA HOUSING AUTH") %>% 
  dplyr::select(-contractorname, -contractoraddress1, -contractoraddress2, -contractorcity, -contractorstate, -contractorzip, -mostrecentinsp) %>% 
  dplyr::mutate(Permit_Year = lubridate::year(permitissuedate), 
                Permit_Month = lubridate::month(permitissuedate)
  ) %>% 
  filter(Permit_Year >= 2010) %>%
  filter(geocode_x != "NA") %>% 
  filter(opa_account_num != "NA") %>% 
  distinct(opa_account_num, permitissuedate, .keep_all = TRUE) %>% 
  st_transform(crs = 3365)

NewCon_list <- NewCon %>% dplyr::select(Permit_Year, opa_account_num) %>% 
  mutate(countNewCon = 1) %>%
  group_split(Permit_Year, keep = TRUE) 

####- This is the problem code, the below st_join turns all my values into NAs
#so I can't do a grid count or split the data by year. If I change the spatial
#join function to "left = FALSE" I get the desired result, but it drops many of the
#rows, is there a way to do the join where I can keep all the rows while also 
#passing through the values. 

#The dropping of the rows is essentially eliminating grids from my original fishnet
#so I can produce a full grid in ggplot/gganimate

### Testing a purr function by taking only one example 
xx <- NewCon_list[[1]]
### Figure out what you need to do with each example
### Mildly convoluted way to join the fishnet ID to points
# then group by year and fishent cell
# then summarise count (could be mean, max, etc) by fishnet
# also summarise year
# then join that count back to fishnet by id
# then turn the NA to a count of zero
# also turn missing year into the correct year

x3 <- st_join(xx, fishnet) %>% 
  group_by(Permit_Year, uniqueID) %>% 
  summarise(cnt = n()) %>% 
  st_drop_geometry() %>% 
  left_join(fishnet, ., by = "uniqueID") %>% 
  mutate(cnt = ifelse(is.na(cnt),0,cnt),
         Permit_Year = ifelse(is.na(Permit_Year),
                            unique(xx$Permit_Year),
                            Permit_Year))

### That works for a single example, but we need it for all
# make a function that can be called from map/apply
point_to_poly <- function(.point, .poly){
  poly_cnt <- st_join(.point, .poly) %>% 
    group_by(Permit_Year, uniqueID) %>% 
    summarise(cnt = n()) %>% 
    st_drop_geometry() %>% 
    left_join(fishnet, ., by = "uniqueID") %>% 
    mutate(cnt = ifelse(is.na(cnt),0,cnt),
           Permit_Year = ifelse(is.na(Permit_Year),
                                unique(.point$Permit_Year),
                                Permit_Year))
}
### test it out on a single example
x4 <- point_to_poly(xx, fishnet)
### make sure it returns the same results
all.equal(x3, x4)

### Try it out in purr for all data
NewConGrid <- map(.x = NewCon_list, ~ point_to_poly(.x, fishnet)) 
### Test to make sure it matches the initial example
all.equal(x3, NewConGrid[[1]])


# NewConGrid <- map(.x = NewCon_list,  ~ st_join(x = fishnet, y = .x, left = TRUE)) 

CombinedConGrid <- mapedit:::combine_list_of_sf(NewConGrid) %>% 
  mutate(Permit_Year = as.integer(Permit_Year))

###CombinedConGrid Output 

PermitGrowth <- ggplot() +
  geom_sf(data=CombinedConGrid, aes(fill = cnt), color = NA) +
  geom_sf(data = neighborhoods, fill = NA, color = "gray65", 
          size = 0.25) +
  scale_fill_viridis_c(begin = 0.2, option = "A",
                       direction = -1, name = "permit\ncount") +
  mapTheme() +
  labs(title = "Count of Construction Permits by Year",
       subtitle = 'Year: {frame_time}') +
  transition_time(Permit_Year) +
  ease_aes('linear')

PermitGrowth