R Notebook

library(sp)
library(spdep)
library(tmap)
library(classInt)
library(grid)
library(gridExtra)
library(lattice)

# data should be clean, no NA no infinite values
# Color scheme ------------------------------------------------------------
bvColors <- c("#e8e8e8","#e4acac","#c85a5a",
           "#b0d5df","#ad9ea5","#985356",
           "#64acbe","#627f8c","#574249")

# loop credit to Stefano De Sabbata
# you can change the colors and the the size of the legend box

# Print bivariate map -----------------------------------------------------

bivariate_choropleth <- function (

    # Function parameters
    bivmap_dataset,         # A SpatialPoligonDataFrame
    bivmap_vars,            # A vector of characters containing the name of the two variables
    bivmap_labels=NA,       # A vector of characters containing the labels for the two variables, to use in the legend
    bivmap_style='quantile',# Classification type for the bins
    bivmap_scale=FALSE      # Use a scale bar

  ) {
  
  bivmap <- get_bivariate_choropleth(
    # Passs parameters on
    # except labels
    bivmap_dataset,
    bivmap_vars,
    bivmap_style,
    bivmap_scale
  )
  
  if (is.na(bivmap_labels)){
    bivmap_labels <- bivmap_vars
  }
  
  # Print map
  suppressWarnings(print( bivmap ))
  
  # Create the square legend
  #x=.9, y=.3, width=.4, height=.4
  # play around with these to figure out which one works the best in your plot
  vp <- viewport(x=.7, y=.12, width=.3, height=.2)
  pushViewport(vp)
  print(levelplot(
    matrix(1:9, nrow=3), 
    axes=FALSE, 
    col.regions=bvColors,
    xlab=list(label=bivmap_labels[1],cex=0.8), 
    ylab=list(label=bivmap_labels[2],cex=0.8), 
    cuts=8, 
    colorkey=FALSE,
    scales=list(draw=0)),
    newpage=FALSE)
  
  # Pop viewport
  popViewport()
}

# Create bivariate map ----------------------------------------------------
# This function actually creates the bivariate map using tmap

get_bivariate_choropleth <- function (
  
  # Function parameters
  bivmap_dataset,         # A SpatialPoligonDataFrame
  bivmap_vars,            # A vector of characters containing the name of the two variables
  bivmap_style='quantile',# Classification type for the bins
  bivmap_scale=FALSE      # Use a scale bar
  
) {
  
  
  # Extract the two specified colums
  # excluding rows with na and infinite values
  
  bivmap_sdf <- bivmap_dataset[, bivmap_vars]
  
  # Renaming the variables to simplify the code below
  colnames(bivmap_sdf@data) <- c("xvar","yvar")
  
  # Create the 3-class categorization per each variable
  bivmap_sdf$xcat <- findCols(classIntervals( bivmap_sdf$xvar, n=3, bivmap_style))
  cat(bivmap_vars[1], "breaks (x-axis):\n")
  print(classIntervals( bivmap_sdf$xvar, n=3, bivmap_style))
  #
  bivmap_sdf$ycat <- findCols(classIntervals( bivmap_sdf$yvar, n=3, bivmap_style))
  cat(bivmap_vars[2], "breaks (y-axis):\n")
  print(classIntervals( bivmap_sdf$yvar, n=3, bivmap_style))
  
  # Combine the above into one 9-class categorization
  bivmap_sdf$bicat <- bivmap_sdf$xcat + (3 * (bivmap_sdf$ycat - 1))
  
  bivmap_sdf$bicol <- bvColors[bivmap_sdf$bicat]
  bivmap_sdf$bicol <- ifelse(is.na(bivmap_sdf$bicol), "#bdbdbd", bivmap_sdf$bicol)
  
  #View(bivmap_sdf@data)
  #View(cbind(bivmap_sdf@data, bivmap_dataset@data))
  
  # Create the map
  bivmap <- tm_shape(bivmap_sdf) + 
    # Fill
    tm_fill(
      "bicol") +
    tm_borders("black", lwd = 0.5) +
    # Remove frame
    tm_layout(frame=FALSE) +
    # Add rhe legend
    tm_legend(scale=0.75)
  
  if (bivmap_scale) {
    bivmap <- bivmap  +
      # Add scale bar
      tm_scale_bar(
        width=0.30,
        position=c("left","bottom"))
  }
  
  # Return bivariate map
  bivmap

}

# so pretty :)

# code
library(sf)
shp <- st_read("/Users/saina/Desktop/311 analysis results/paper/First paper/floodANDNewpermit.shp")

## Reading layer `floodANDNewpermit' from data source 
##   `/Users/saina/Desktop/311 analysis results/paper/First paper/floodANDNewpermit.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 305 features and 15 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -71.29659 ymin: 42.22793 xmax: -70.95839 ymax: 42.46563
## Geodetic CRS:  NAD83

world_spdf <- as(shp, "Spatial")
tmap_mode("plot")
# Plot bivariate choroplet map
bivariate_choropleth(world_spdf, c("Fld_Rsk", "nw_cns_"), bivmap_labels = c("risk\n low  >>  high", "New const.\n low  >>  high"))

## Fld_Rsk breaks (x-axis):
## style: quantile
##   one of 903 possible partitions of this variable into 3 classes
##  [0,0)  [0,6) [6,48] 
##      0    190    115 
## nw_cns_ breaks (y-axis):
## style: quantile
##          [0,0.0382) [0.0382,0.07878333) [0.07878333,0.2857] 
##                 101                 102                 102

#world_spdf <- as(shp, "Spatial")

# Plot bi variate choroplet map
#bivariate_choropleth(world_spdf, c("Flood_Risk", "nw_cnst_new"), bivmap_labels = c("risk\n low  >>  high", "alterations\n low  >>  high"))