sp objects

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

suppressWarnings(library(SfSpHelpers))
suppressWarnings(library(raster))
suppressWarnings(library(sf))
suppressWarnings(library(magrittr))
suppressWarnings(library(dplyr))
suppressWarnings(library(ggplot2))
suppressWarnings(library(RColorBrewer))

Data

  #Trees
  shp_trees <-  st_read('https://www.donneesquebec.ca/recherche/dataset/bc5afddf-9439-4e96-84fb-f91847b722be/resource/bbdca0dd-82df-42f9-845b-32348debf8ab/download/vdq-arbrepotentielremarquable.geojson')

  #Neighborhoods
  shp_neigh <-  get_zipped_remote_shapefile("https://www.donneesquebec.ca/recherche/dataset/5b1ae6f2-6719-46df-bd2f-e57a7034c917/resource/508594dc-b090-407c-9489-73a1b46a8477/download/vdq-quartier.zip")

2D raster KDE plot - default

  rasterKDECentroids <- st_kde(shp_trees , cellsize = 0.001, bandwith =c(.001, .001 ) )
  plot(rasterKDECentroids$layer, main='Remarkable tree density - Quebec City')

GGplot version

  df_rast <- rasterKDECentroids$layer %>% as.data.frame(xy=T)

  if(! ('NOM' %in% colnames(shp_neigh))){
    knitr::knit_exit()
  }
  extended_cols <- colorRampPalette(brewer.pal(8, "Set3"))(n_distinct(shp_neigh[['NOM']]))

  ggplot() + 
    geom_tile(data=df_rast, aes(x=x,y=y,fill=layer),alpha=0.2) + 
    geom_sf(data=shp_neigh, aes(col=NOM),alpha=0) + 
    scale_fill_viridis_c() + 
    scale_color_manual(values = extended_cols) + 
    ggplot2::theme_dark(base_family="Roboto Condensed", base_size=11.5) +
    guides( color = "none")

Zoom in using functions in bbox

shp_bbox_raster <- bbox_from_vector( extent(rasterKDECentroids) , crs=4326)

shp_neigh_filtered <- st_intersection(shp_neigh %>% st_zm(),
                                      shp_bbox_raster)

  ggplot() + 
    geom_tile(data=df_rast, aes(x=x,y=y,fill=layer),alpha=0.2) + 
    geom_sf(data=shp_neigh_filtered, aes(col=NOM),alpha=0) + 
    scale_fill_viridis_c() + 
    scale_color_manual(values = extended_cols) + 
    ggplot2::theme_dark(base_family="Roboto Condensed", base_size=11.5) +
    guides( color = 'none')

Polygon contours

  #get_polygon_heatmap Works with polygons also (takes centroid implicitely)
  shp_polyons <- get_polygon_heatmap(shp_trees , bw=.001, gsize=500 )


#Can use the colors produced automatically, but this is a red to yellow gradient 
ggplot(shp_polyons)+
    geom_sf(aes( fill = colors),lwd=0) +
    scale_fill_identity() +
    geom_sf(data=shp_neigh, aes(col=NOM),alpha=0) + 
    ggplot2::theme_minimal(base_family="Roboto Condensed", base_size=11.5) +
    guides( color = 'none') 

#Can simply use viridis discrete
ggplot(shp_polyons)+
  geom_sf( aes(fill=colors) ,lwd=0) +
  geom_sf(data=shp_neigh, aes(col=NOM),alpha=0) + 
  scale_fill_viridis_d()+
  coord_sf(datum = NA) + 
  ggplot2::theme_minimal(base_family="Roboto Condensed", base_size=11.5) +
  guides( color = 'none', fill = 'none') +
  ggtitle("Quebec city remakable tree density") + 
  labs(subtitle = "2D kernel density estimate",
         caption ='Source: https://www.donneesquebec.ca - vdq-arbrepotentielremarquable.geojson')

ggsave('tree_kde.png', width = 7, height = 7)