Note:

The Map shows the major municipalities in British Columbia, Canada. The population and area of the cities are shown in blue and red color circles, respectively (need to zoom). The size of the circles represent the relative size of the population and area of any municipality. The popup marker also display the name, population and area of the city.

Appendix (Code):

Load required libraries.

library(htmltab)
library(dplyr)
library(leaflet)
## Read Table from wikipedia

table1 <- htmltab(doc = "https://en.wikipedia.org/wiki/List_of_municipalities_in_British_Columbia#Municipalities", 1)

## Take a subset of the data table and rename the columns (contains city, population, area)

BC <- subset(table1, 
             select = c("Name", 
                        "Status", 
                        "2016 Census of Population >> Population (2016)", 
                        "2016 Census of Population >> Land area (km²)"))
BC <- rename(BC, Population="2016 Census of Population >> Population (2016)")
BC <- rename(BC, Area="2016 Census of Population >> Land area (km²)")

## Cleaning the data

BC<- BC[-c(162:172),]
BC$Population <- as.numeric(gsub(",", "", BC$Population))
BC$Area <- as.numeric(gsub(",", "", BC$Area))
BC <- BC[-which(is.na(BC$Population)),]

## Read the 2016 census report and take a subset  (contains city, lat and lng)

bccor <- read.csv("2016_92-151_XBB.csv")
BC1 <- subset(bccor, PRename.PRanom == "British Columbia", 
                 select = c("CSDname.SDRnom", 
                            "DArplat.Adlat", 
                            "DArplong.ADlong"))

## match with the previous data frame and clean the data

match1 <- BC1$CSDname.SDRnom %in% BC$Name
df <- BC1[match1,]
df <- df[, c("CSDname.SDRnom", "DArplat.Adlat", "DArplong.ADlong")]
df <- df[!duplicated(df[1]),]
colnames(df)[1] <- "Name"

## Merge two data frame based on city name.

df <- merge(BC, df, by.x = "Name", all.x = TRUE)
df <- df[!duplicated(df[1]),]
df <- df[!is.na(df$DArplat.Adlat),]

colnames(df)[5] <- "latitude"
colnames(df)[6] <- "longitude"

df_latlong <- df[, c("latitude", "longitude")]
## The final code chunk for map

df_latlong %>% 
  leaflet() %>%
  addTiles() %>%
  addMarkers(clusterOptions = markerClusterOptions(maxClusterRadius = 15), 
             popup = paste("City:", df$Name, "<br>",
                           "Population:", df$Population, "<br>",
                           "Area (km^2):", df$Area)) %>%
  addCircles(weight = 1, radius = sqrt(df$Population) * 5, color = "blue") %>%
  addCircles(weight = 1, radius = sqrt(df$Area) * 100, color = "red") %>%
  addLegend(labels = c("Population", "Area"), colors = c("blue", "red"))