SgPopulation
Daniel
8/25/2020
packages = c('rgdal', 'spdep', 'ClustGeo', 'tmap', 'sf', 'ggpubr', 'cluster', 'heatmaply', 'corrplot', 'psych', 'tidyverse', "factoextra", "NbClust", "psych")
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}## Loading required package: rgdal## Loading required package: sp## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-2## Loading required package: spdep## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: sf## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0## Loading required package: ClustGeo## Loading required package: tmap## Loading required package: ggpubr## Loading required package: ggplot2## Loading required package: cluster## Loading required package: heatmaply## Loading required package: plotly##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layout## Loading required package: viridis## Loading required package: viridisLite## Registered S3 method overwritten by 'seriation':
## method from
## reorder.hclust gclus##
## ======================
## Welcome to heatmaply version 1.1.0
##
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
##
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## Or contact: <tal.galili@gmail.com>
## ======================## Loading required package: corrplot## corrplot 0.84 loaded## Loading required package: psych##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## Loading required package: tidyverse## ── Attaching packages ─────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ tibble 2.1.3 ✓ dplyr 0.8.4
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ✓ purrr 0.3.3## ── Conflicts ────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x psych::%+%() masks ggplot2::%+%()
## x psych::alpha() masks ggplot2::alpha()
## x dplyr::filter() masks plotly::filter(), stats::filter()
## x dplyr::lag() masks stats::lag()## Loading required package: factoextra## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa## Loading required package: NbClustImport Geospatial Layer
mpsz <- st_read(dsn = "dataSource/geospatial",
layer = "MP14_SUBZONE_WEB_PL")## Reading layer `MP14_SUBZONE_WEB_PL' from data source `/Users/daniel/Google Drive (daniel.soh.2017@smu.edu.sg)/SMU | GD/1. MGMT334 Social Entrepreneurship Practicum/Group/Data/dataSource/geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## proj4string: +proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defsmpsz_invalid <- mpsz %>% filter(!st_is_valid(mpsz))
st_is_valid(mpsz_invalid, reason = TRUE)## [1] "Ring Self-intersection[27932.3925999999 21982.7971999999]"
## [2] "Ring Self-intersection[26885.4439000003 26668.3121000007]"
## [3] "Ring Self-intersection[26920.1689999998 26978.5440999996]"
## [4] "Ring Self-intersection[15432.4749999996 31319.716]"
## [5] "Ring Self-intersection[12861.3828999996 32207.4923]"
## [6] "Ring Self-intersection[19681.2353999997 31294.4521999992]"
## [7] "Ring Self-intersection[41375.108 40432.8588999994]"
## [8] "Ring Self-intersection[38542.2260999996 44605.4089000002]"
## [9] "Ring Self-intersection[21702.5623000003 48125.1154999994]"Remove invalid
mpsz_valid <- st_make_valid(mpsz) %>% filter(!st_is_valid(mpsz))Check valid
mpsz <- st_make_valid(mpsz)
any(st_is_valid(mpsz))## [1] TRUEConvert CRS to WGS84
mpsz <- st_transform(mpsz, crs=4326)merged_sf <- mpsz %>%
mutate(area = st_area(mpsz)) %>%
arrange(SUBZONE_N) %>%
select(SUBZONE_N, area)population <-read_csv("dataSource/aspatial/planning-area-subzone-age-group-sex-and-type-of-dwelling-june-2011-2019.csv")## Parsed with column specification:
## cols(
## planning_area = col_character(),
## subzone = col_character(),
## age_group = col_character(),
## sex = col_character(),
## type_of_dwelling = col_character(),
## resident_count = col_double(),
## year = col_double()
## )Filter population data to 2019 and sum of resident_count by dwelling in each subzone Only the latest population data(2019) is used. The sum of resident_count in all sex and age_group is used to calculate the population by dwelling. From this will we get the following:
HDB1-2RM dwellers HDB3-4RM dwellers HDB5RM-Ec dweller Condominium and apartment dwellers Landed property dwellers ###
population_2019_dwelling <- population %>%
# Filter year 2019
filter(year==2019)%>%
select(2:6) %>%
# Merge sex
pivot_wider(names_from=c(sex), values_from = resident_count,values_fill=list(resident_count=0)) %>%
mutate(population = rowSums(.[4:5])) %>%
select(1:3,6) %>%
# Merge age_group
pivot_wider(names_from=c(age_group), values_from = population,values_fill=list(population=0)) %>%
mutate(population = rowSums(.[3:21])) %>%
select(1:2,22) %>%
# Pivot Type of Dwelling to Columns
pivot_wider(names_from=c(type_of_dwelling), values_from = population,values_fill=list(population=0)) %>%
select(1:5,7:8) %>%
# Merge HDB 3-Room Flats and HDB 4-Room Flats
mutate("HDB 3- and 4-Room Flats" = rowSums(.[3:4])) %>%
select(1:2,8,5:7) %>%
# Make subzone Uppercase
mutate(subzone = toupper(subzone))population_2019_age <- population %>%
# Filter year 2019
filter(year==2019)%>%
select(2:6) %>%
# Merge sex
pivot_wider(names_from=c(sex), values_from = resident_count,values_fill=list(resident_count=0)) %>%
mutate(population = rowSums(.[4:5])) %>%
select(1:3,6) %>%
# Merge Type of Dwelling
pivot_wider(names_from=c(type_of_dwelling), values_from = population,values_fill=list(population=0))%>%
mutate(population = rowSums(.[3:10])) %>%
select(1:2,11) %>%
# Pivot Age Group to Columns
pivot_wider(names_from=c(age_group), values_from = population,values_fill=list(population=0)) %>%
# Calculate population in selected categories
mutate(economy_active = rowSums(.[7:14])) %>%
mutate(young_group = rowSums(.[2:6])) %>%
mutate(aged_group = rowSums(.[15:20])) %>%
mutate(population = rowSums(.[2:20])) %>%
select(1,21:24) %>%
# Make subzone Uppercase
mutate(subzone = toupper(subzone))merged_sf <- left_join(merged_sf, population_2019_age, by=c("SUBZONE_N"="subzone"))## Warning: Column `SUBZONE_N`/`subzone` joining factor and character vector,
## coercing into character vectormerged_sf <- left_join(merged_sf, population_2019_dwelling, by=c("SUBZONE_N"="subzone"))
merged_sf## Simple feature collection with 323 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
## CRS: EPSG:4326
## First 10 features:
## SUBZONE_N area economy_active young_group aged_group
## 1 ADMIRALTY 1261648.7 [m^2] 8380 4370 1360
## 2 AIRPORT ROAD 494503.4 [m^2] 0 0 0
## 3 ALEXANDRA HILL 1030378.6 [m^2] 7560 2910 3310
## 4 ALEXANDRA NORTH 293706.4 [m^2] 1420 580 120
## 5 ALJUNIED 2959365.3 [m^2] 24330 8350 7510
## 6 ANAK BUKIT 2768350.3 [m^2] 12490 5910 3850
## 7 ANCHORVALE 1499107.2 [m^2] 27740 15250 3620
## 8 ANG MO KIO TOWN CENTRE 316882.0 [m^2] 2790 1360 740
## 9 ANSON 103238.5 [m^2] 0 0 0
## 10 BALESTIER 1926621.8 [m^2] 19320 7350 6090
## population HDB 1- and 2-Room Flats HDB 3- and 4-Room Flats
## 1 14110 1140 6750
## 2 0 0 0
## 3 13780 3910 6320
## 4 2120 0 0
## 5 40190 2120 18020
## 6 22250 160 2090
## 7 46610 1680 17280
## 8 4890 0 1120
## 9 0 0 0
## 10 32760 1790 15490
## HDB 5-Room and Executive Flats Landed Properties
## 1 6220 0
## 2 0 0
## 3 3120 0
## 4 0 0
## 5 6400 460
## 6 3170 7600
## 7 22650 0
## 8 1830 0
## 9 0 0
## 10 4660 570
## Condominiums and Other Apartments geometry
## 1 0 MULTIPOLYGON (((103.8285 1....
## 2 0 MULTIPOLYGON (((103.9014 1....
## 3 0 MULTIPOLYGON (((103.8144 1....
## 4 2120 MULTIPOLYGON (((103.8174 1....
## 5 11930 MULTIPOLYGON (((103.8913 1....
## 6 9020 MULTIPOLYGON (((103.771 1.3...
## 7 5000 MULTIPOLYGON (((103.896 1.3...
## 8 1930 MULTIPOLYGON (((103.8485 1....
## 9 0 MULTIPOLYGON (((103.8441 1....
## 10 9610 MULTIPOLYGON (((103.8623 1....plotMap <- function(sf_df, sf_str) {
tm_shape(sf_df)+
tm_layout(legend.outside = TRUE, legend.hist.width = 2)+
tm_polygons(col=sf_str, legend.hist = TRUE)
}population_map <- plotMap(merged_sf, "population")
young_map <- plotMap(merged_sf, "young_group")
economy_active_map <- plotMap(merged_sf, "economy_active")
aged_group_map <- plotMap(merged_sf, "aged_group")
tmap_arrange(population_map, young_map, economy_active_map, aged_group_map)## Some legend labels were too wide. These labels have been resized to 0.51, 0.51, 0.51, 0.47, 0.44, 0.44. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Some legend labels were too wide. These labels have been resized to 0.55, 0.51, 0.51, 0.51, 0.51, 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Some legend labels were too wide. These labels have been resized to 0.51, 0.51, 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Some legend labels were too wide. These labels have been resized to 0.55, 0.51, 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
HDB12_map <- plotMap(merged_sf, "HDB 1- and 2-Room Flats")
HDB34_map <- plotMap(merged_sf, "HDB 3- and 4-Room Flats")
HDB5Ec_map <- plotMap(merged_sf, "HDB 5-Room and Executive Flats")
CondoApt_map <- plotMap(merged_sf, "Condominiums and Other Apartments")
Landed_map <- plotMap(merged_sf, "Landed Properties")
tmap_arrange(HDB12_map, HDB34_map, HDB5Ec_map)## Some legend labels were too wide. These labels have been resized to 0.59, 0.59, 0.59, 0.59. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Some legend labels were too wide. These labels have been resized to 0.51, 0.51, 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Some legend labels were too wide. These labels have been resized to 0.51, 0.51, 0.51, 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tmap_arrange(CondoApt_map, Landed_map, ncol=1)
tmap_mode("view")## tmap mode set to interactive viewingtm_shape(merged_sf) +
tm_polygons("HDB 1- and 2-Room Flats")tm_shape(merged_sf) +
tm_polygons("HDB 1- and 2-Room Flats")tm_shape(merged_sf) +
tm_polygons("HDB 3- and 4-Room Flats")tm_shape(merged_sf) +
tm_polygons("HDB 5-Room and Executive Flats")tm_shape(merged_sf) +
tm_polygons("Condominiums and Other Apartments")