R/nlchoropleth.R

In MHaringa/nlmaps: Choropleth maps

#' Create choropleth simple feature object from points
#'
#' @description Create choropleth simple feature object from points.
#'
#' @param sf_map Simple features object. Note that shapefiles consist of more than one file, all with identical basename, which should reside in the same directory.
#' @param df Data.frame of interest. The data.frame should contain the column names 'lon' and 'lat'
#' @param oper An operation on the polygon level.
#' @param crs Coordinate reference system: integer with the EPSG code, or character with proj4string
#'
#' @return standard feature object
#' @export choropleth_sf
#'
#' @import sf
#' @import dplyr
#'
#' @author Martin Haringa
#'
#' @examples
#' choropleth_sf(nl_postcode2, insurance, sum(amount, na.rm = TRUE))
#' \dontrun{
#' shp_read <- sf::st_read(~/path/to/file.shp)
#' choropleth_sf(shp_read, insurance, sum(amount, na.rm = TRUE))
#'}
choropleth_sf <- function(sf_map, df, oper, crs = 4326){

  oper <- enquo(oper)

  # https://www.r-spatial.org/r/2018/10/25/ggplot2-sf.html
  shp_wgs84 <- sf_map %>%
    sf::st_buffer(0) %>% # Make invalid geometries valid
    st_transform(crs = crs) %>% # Convert coordinates to WGS84
    dplyr::mutate(id = 1:nrow(sf_map))

  if( !all(c("lon", "lat") %in% names(df)) ) {stop("Data.frame should contain column names 'lon' and 'lat'.")}

  df_sf <- sf::st_as_sf(df, coords = c("lon", "lat"), crs = crs)

  suppressMessages({
    df_map_sf <- sf::st_join(shp_wgs84, df_sf)
  })

  # Change from sf to data.frame
  df_map <- df_map_sf
  sf::st_geometry(df_map) <- NULL

  # Aggregate
  df_map_sf2 <- df_map %>%
    dplyr::group_by(id) %>%
    dplyr::summarize(output = !! oper) %>%
    dplyr::ungroup()

  out <- dplyr::left_join(shp_wgs84, df_map_sf2, by = "id")
  return(out)
}