llibrary(sf)
library(tidyverse)
rm(list=ls())

# load raster of NLCD with terra
nlcd <- terra::rast('Documents/GIS.Data/agriculture/nlcd2019_ny.tif')

# load counties in sf using tigris
cos <- tigris::counties('ny')

# exactextractr is a fast way to do zonal statistics and can be used
# for zonal histograms using this code
#
# append_cols = T  ## include columns from shapefile to allow rejoining at end
# summarize_df = T ## pass the dataframecreated by extact_extract to the function
# function = .. calculate the fraction of the area using frac (fraction of area) option in exactexractr

nlcd_co <- exactextractr::exact_extract(nlcd, cos,
                                        append_cols=T,
                                        fun='frac' )
                                    
# this next step is optional, it just takes  long format and pivots_wider
# for an output similiar to what you would get in ArcGIS or QGIS
# then inner join restores spatial data by joining on matching rows
# also we convert to acres, which is more useful then square meters.

nlcd_co %>%
    mutate(area = units::set_units(area,'acres')) %>%
    pivot_wider(names_from = 'nlcd_code', values_from = 'area',
        names_glue = "nlcd{nlcd_code}" ) %>%
    inner_join(cos, .) -> nlcd_co

# Create the Table Shown Below
nlcd_co %>% 
  arrange(NAME) %>%
  select(NAME, starts_with('nlcd'), -nlcdNaN) %>%
  mutate(across(starts_with('nlcd'), ~round(.x))) %>%
  st_drop_geometry()