IS415_Take_Home_Exercise_alt
leeSunho
9/21/2020
Rpub link
Objective
Analyse the geographic distribution childcare services and compare the chilcare between 2017 to 2020. Use appropriate map analysis and spatial point patterns analysis technique.
Understand of the Singpaore Zone
Singapore regions are sub divided into 55 planning area and this includes two water catchment areas. the largest area is the West Region and Central Region is the most populous area.
knitr::include_graphics("picture/sg.png")
Data Preperation, Importing and Wrangling
Before I begin the analysis, I will be spending time preparing datasets for the requirements.
importing the required packages
Firstly, I imported necessary packages. R packages are a collection of R functions which contain code and some sample data. and They are stores under a library in R server. Different library have different purposes.
packages = c('sf', 'tmap', 'tidyverse','dplyr', 'magrittr','rgdal','maptools', 'raster', 'spatstat','rgdal')
for (p in packages){
if(!require(p, character.only = T)){
install.packages(p)
}
library(p,character.only = T)
}Import datasets
Master Plan 2014 Planning Area Boundary
the dataset is imported from: https://data.gov.sg/dataset/master-plan-2014-subzone-boundary-web
mpsz <- st_read(dsn = "data/2017/geospatial",
layer = "MP14_SUBZONE_WEB_PL")## Reading layer `MP14_SUBZONE_WEB_PL' from data source `C:\Users\sunho\Documents\SMU\3.2\gis\Take home\Take_Home\data\2017\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
## projected CRS: SVY21Examine the contents of mpsz
Check for CRS
st_crs(mpsz)## Coordinate Reference System:
## User input: SVY21
## wkt:
## PROJCRS["SVY21",
## BASEGEOGCRS["SVY21[WGS84]",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]As mpsz epsg is 9001 convert into 3414
mpsz <- st_set_crs(mpsz, 3414)Check again for mpsz
st_crs(mpsz)## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]Check for NA data
mpsz[rowSums(is.na(mpsz))!=0,]## Simple feature collection with 0 features and 15 fields
## bbox: xmin: NA ymin: NA xmax: NA ymax: NA
## projected CRS: SVY21 / Singapore TM
## [1] OBJECTID SUBZONE_NO SUBZONE_N SUBZONE_C CA_IND PLN_AREA_N
## [7] PLN_AREA_C REGION_N REGION_C INC_CRC FMEL_UPD_D X_ADDR
## [13] Y_ADDR SHAPE_Leng SHAPE_Area geometry
## <0 rows> (or 0-length row.names)Check validity of geometrics
mpsz_check<- st_is_valid(mpsz)
length(which(mpsz_check == FALSE))## [1] 9since there are 9 invalid geometries, make it into valid
mpsz<- st_make_valid(mpsz)
mpsz_check <- st_is_valid(mpsz)
length(which(mpsz_check == FALSE))## [1] 0Plot the spatial data and understand different subzone
plot(mpsz)
Import aspatial data
Importing aspatial dataet which includes Planning AreaSubzone, Age Group, Sex and Type of Dwelling
Import the aspatial data “Singapore Residents by Planning AreaSubzone, Age Group, Sex and Type of Dwelling” Dataset can be found: https://www.singstat.gov.sg/find-data/search-by-theme/population/geographic-distribution/latest-data
popdata <- read_csv("data/2017/aspatial/respopagesextod2011to2019.csv")Check the population dataset to ensure it is imported correctly
list(popdata)## [[1]]
## # A tibble: 883,728 x 7
## PA SZ AG Sex TOD Pop Time
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males HDB 1- and 2-Room Flats 0 2011
## 2 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males HDB 3-Room Flats 10 2011
## 3 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males HDB 4-Room Flats 30 2011
## 4 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males HDB 5-Room and Executive~ 50 2011
## 5 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males HUDC Flats (excluding th~ 0 2011
## 6 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males Landed Properties 0 2011
## 7 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males Condominiums and Other A~ 40 2011
## 8 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Males Others 0 2011
## 9 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Femal~ HDB 1- and 2-Room Flats 0 2011
## 10 Ang Mo K~ Ang Mo Kio Town~ 0_to~ Femal~ HDB 3-Room Flats 10 2011
## # ... with 883,718 more rowsData Wrangling
Child Care Centers (CCCs) provide full day and half-day care programs to children below the age of 7 years. Some centers provide infant care services to children aged 2 months to 18 months. Some age-related development periods and examples of defined intervals include: newborn (ages 0–4 weeks); infant (ages 4 weeks – 1 year); toddler (ages 12 months-24 months); preschooler (ages 2–5 years); school-aged child (ages 6–12 years); adolescent (ages 13–19).
However, the dataset I have is grouped by up to 0 to 4 and 4 to 10. Since, Hence, I will be using 0_to4 age group.
I have filtered the year to 2019 and 2017 separately for further analysis later. As well as, made a special row for Child.
2019
mpszpop2019 <- popdata %>%
filter(Time == 2019) %>%
group_by(PA,SZ,AG) %>%
summarize(`Pop` = sum(`Pop`)) %>%
ungroup() %>%
spread(AG, Pop) %>%
mutate(Child = `0_to_4`) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[9:13])+
rowSums(.[15:17]))%>%
mutate(`AGED`=rowSums(.[18:22])) %>%
mutate(`TOTAL`=rowSums(.[5:22])) %>%
mutate(`DEPENDENCY` = (`Child` + `AGED`)
/`ECONOMY ACTIVE`) %>%
mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper)) %>%
dplyr::select(`PA`, `SZ`, `Child`, `ECONOMY ACTIVE`, `AGED`,
`TOTAL`, `DEPENDENCY`) %>%
filter(`ECONOMY ACTIVE` > 0)
mpszpop2019## # A tibble: 234 x 7
## PA SZ Child `ECONOMY ACTIVE` AGED TOTAL DEPENDENCY
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ANG MO KIO ANG MO KIO TOWN CEN~ 170 2600 440 4560 0.235
## 2 ANG MO KIO CHENG SAN 1050 15480 3110 27160 0.269
## 3 ANG MO KIO CHONG BOON 840 14330 3370 25910 0.294
## 4 ANG MO KIO KEBUN BAHRU 730 12520 2690 22120 0.273
## 5 ANG MO KIO SEMBAWANG HILLS 200 3250 740 6360 0.289
## 6 ANG MO KIO SHANGRI-LA 540 8410 1970 15440 0.298
## 7 ANG MO KIO TAGORE 220 3880 790 7610 0.260
## 8 ANG MO KIO TOWNSVILLE 750 11660 2830 20620 0.307
## 9 ANG MO KIO YIO CHU KANG EAST 160 2060 410 3970 0.277
## 10 ANG MO KIO YIO CHU KANG WEST 760 13220 2220 23300 0.225
## # ... with 224 more rows2017
mpszpop2017 <- popdata %>%
filter(Time == 2017) %>%
group_by(PA,SZ,AG) %>%
summarize(`Pop` = sum(`Pop`)) %>%
ungroup() %>%
spread(AG, Pop) %>%
mutate(Child = `0_to_4`) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[9:13])+
rowSums(.[15:17]))%>%
mutate(`AGED`=rowSums(.[18:22])) %>%
mutate(`TOTAL`=rowSums(.[5:22])) %>%
mutate(`DEPENDENCY` = (`Child` + `AGED`)
/`ECONOMY ACTIVE`) %>%
mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper)) %>%
dplyr::select(`PA`, `SZ`, `Child`, `ECONOMY ACTIVE`, `AGED`,
`TOTAL`, `DEPENDENCY`) %>%
filter(`ECONOMY ACTIVE` > 0)
mpszpop2017## # A tibble: 233 x 7
## PA SZ Child `ECONOMY ACTIVE` AGED TOTAL DEPENDENCY
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ANG MO KIO ANG MO KIO TOWN CEN~ 190 2690 450 4640 0.238
## 2 ANG MO KIO CHENG SAN 1120 15590 2970 27250 0.262
## 3 ANG MO KIO CHONG BOON 850 14540 2970 25870 0.263
## 4 ANG MO KIO KEBUN BAHRU 680 12190 2480 21690 0.259
## 5 ANG MO KIO SEMBAWANG HILLS 190 3290 690 6400 0.267
## 6 ANG MO KIO SHANGRI-LA 590 8940 1920 16150 0.281
## 7 ANG MO KIO TAGORE 280 3920 790 7720 0.273
## 8 ANG MO KIO TOWNSVILLE 960 12220 2910 21750 0.317
## 9 ANG MO KIO YIO CHU KANG EAST 170 1970 430 3880 0.305
## 10 ANG MO KIO YIO CHU KANG WEST 920 13750 2290 24390 0.233
## # ... with 223 more rowsCheck filtered dataframe for 2019
list(mpszpop2019)## [[1]]
## # A tibble: 234 x 7
## PA SZ Child `ECONOMY ACTIVE` AGED TOTAL DEPENDENCY
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ANG MO KIO ANG MO KIO TOWN CEN~ 170 2600 440 4560 0.235
## 2 ANG MO KIO CHENG SAN 1050 15480 3110 27160 0.269
## 3 ANG MO KIO CHONG BOON 840 14330 3370 25910 0.294
## 4 ANG MO KIO KEBUN BAHRU 730 12520 2690 22120 0.273
## 5 ANG MO KIO SEMBAWANG HILLS 200 3250 740 6360 0.289
## 6 ANG MO KIO SHANGRI-LA 540 8410 1970 15440 0.298
## 7 ANG MO KIO TAGORE 220 3880 790 7610 0.260
## 8 ANG MO KIO TOWNSVILLE 750 11660 2830 20620 0.307
## 9 ANG MO KIO YIO CHU KANG EAST 160 2060 410 3970 0.277
## 10 ANG MO KIO YIO CHU KANG WEST 760 13220 2220 23300 0.225
## # ... with 224 more rowsCheck filtered dataframe for 2017
list(mpszpop2017)## [[1]]
## # A tibble: 233 x 7
## PA SZ Child `ECONOMY ACTIVE` AGED TOTAL DEPENDENCY
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ANG MO KIO ANG MO KIO TOWN CEN~ 190 2690 450 4640 0.238
## 2 ANG MO KIO CHENG SAN 1120 15590 2970 27250 0.262
## 3 ANG MO KIO CHONG BOON 850 14540 2970 25870 0.263
## 4 ANG MO KIO KEBUN BAHRU 680 12190 2480 21690 0.259
## 5 ANG MO KIO SEMBAWANG HILLS 190 3290 690 6400 0.267
## 6 ANG MO KIO SHANGRI-LA 590 8940 1920 16150 0.281
## 7 ANG MO KIO TAGORE 280 3920 790 7720 0.273
## 8 ANG MO KIO TOWNSVILLE 960 12220 2910 21750 0.317
## 9 ANG MO KIO YIO CHU KANG EAST 170 1970 430 3880 0.305
## 10 ANG MO KIO YIO CHU KANG WEST 920 13750 2290 24390 0.233
## # ... with 223 more rowsFrom the distribution of number of child, we can notice that 2019 have lesser child number compare to 2017. This means that more couples are giving birth.
Georelation Join
Joining mpsz and popdata
let’s join the child population data and singapore map to understand the distribution of child born in Singapore.
2019
mpsz2019 <- left_join(mpsz, mpszpop2019,
by = c("SUBZONE_N" = "SZ"))
mpsz2019## Simple feature collection with 323 features and 21 fields
## geometry type: GEOMETRY
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS: SVY21 / Singapore TM
## First 10 features:
## OBJECTID SUBZONE_NO SUBZONE_N SUBZONE_C CA_IND PLN_AREA_N
## 1 1 1 MARINA SOUTH MSSZ01 Y MARINA SOUTH
## 2 2 1 PEARL'S HILL OTSZ01 Y OUTRAM
## 3 3 3 BOAT QUAY SRSZ03 Y SINGAPORE RIVER
## 4 4 8 HENDERSON HILL BMSZ08 N BUKIT MERAH
## 5 5 3 REDHILL BMSZ03 N BUKIT MERAH
## 6 6 7 ALEXANDRA HILL BMSZ07 N BUKIT MERAH
## 7 7 9 BUKIT HO SWEE BMSZ09 N BUKIT MERAH
## 8 8 2 CLARKE QUAY SRSZ02 Y SINGAPORE RIVER
## 9 9 13 PASIR PANJANG 1 QTSZ13 N QUEENSTOWN
## 10 10 7 QUEENSWAY QTSZ07 N QUEENSTOWN
## PLN_AREA_C REGION_N REGION_C INC_CRC FMEL_UPD_D X_ADDR
## 1 MS CENTRAL REGION CR 5ED7EB253F99252E 2014-12-05 31595.84
## 2 OT CENTRAL REGION CR 8C7149B9EB32EEFC 2014-12-05 28679.06
## 3 SR CENTRAL REGION CR C35FEFF02B13E0E5 2014-12-05 29654.96
## 4 BM CENTRAL REGION CR 3775D82C5DDBEFBD 2014-12-05 26782.83
## 5 BM CENTRAL REGION CR 85D9ABEF0A40678F 2014-12-05 26201.96
## 6 BM CENTRAL REGION CR 9D286521EF5E3B59 2014-12-05 25358.82
## 7 BM CENTRAL REGION CR 7839A8577144EFE2 2014-12-05 27680.06
## 8 SR CENTRAL REGION CR 48661DC0FBA09F7A 2014-12-05 29253.21
## 9 QT CENTRAL REGION CR 1F721290C421BFAB 2014-12-05 22077.34
## 10 QT CENTRAL REGION CR 3580D2AFFBEE914C 2014-12-05 24168.31
## Y_ADDR SHAPE_Leng SHAPE_Area PA Child ECONOMY ACTIVE AGED
## 1 29220.19 5267.381 1630379.3 <NA> NA NA NA
## 2 29782.05 3506.107 559816.2 OUTRAM 170 3580 1130
## 3 29974.66 1740.926 160807.5 SINGAPORE RIVER 0 50 10
## 4 29933.77 3313.625 595428.9 BUKIT MERAH 440 6880 1910
## 5 30005.70 2825.594 387429.4 BUKIT MERAH 520 5810 1270
## 6 29991.38 4428.913 1030378.8 BUKIT MERAH 380 7010 1770
## 7 30230.86 3275.312 551732.0 BUKIT MERAH 520 7970 2110
## 8 30222.86 2208.619 290184.7 SINGAPORE RIVER 0 30 10
## 9 29893.78 6571.323 1084792.3 QUEENSTOWN 170 2380 350
## 10 30104.18 3454.239 631644.3 QUEENSTOWN 0 180 10
## TOTAL DEPENDENCY geometry
## 1 NA NA MULTIPOLYGON (((31495.56 30...
## 2 6410 0.36312849 MULTIPOLYGON (((29092.28 30...
## 3 70 0.20000000 MULTIPOLYGON (((29932.33 29...
## 4 12870 0.34156977 MULTIPOLYGON (((27131.28 30...
## 5 10120 0.30808950 MULTIPOLYGON (((26451.03 30...
## 6 13190 0.30670471 MULTIPOLYGON (((25899.7 297...
## 7 14240 0.32998745 MULTIPOLYGON (((27746.95 30...
## 8 60 0.33333333 MULTIPOLYGON (((29351.26 29...
## 9 4190 0.21848739 MULTIPOLYGON (((20996.49 30...
## 10 280 0.05555556 MULTIPOLYGON (((24472.11 29...2017
mpsz2017 <- left_join(mpsz, mpszpop2017, by = c("SUBZONE_N" = "SZ"))
mpsz2017## Simple feature collection with 323 features and 21 fields
## geometry type: GEOMETRY
## dimension: XY
## bbox: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS: SVY21 / Singapore TM
## First 10 features:
## OBJECTID SUBZONE_NO SUBZONE_N SUBZONE_C CA_IND PLN_AREA_N
## 1 1 1 MARINA SOUTH MSSZ01 Y MARINA SOUTH
## 2 2 1 PEARL'S HILL OTSZ01 Y OUTRAM
## 3 3 3 BOAT QUAY SRSZ03 Y SINGAPORE RIVER
## 4 4 8 HENDERSON HILL BMSZ08 N BUKIT MERAH
## 5 5 3 REDHILL BMSZ03 N BUKIT MERAH
## 6 6 7 ALEXANDRA HILL BMSZ07 N BUKIT MERAH
## 7 7 9 BUKIT HO SWEE BMSZ09 N BUKIT MERAH
## 8 8 2 CLARKE QUAY SRSZ02 Y SINGAPORE RIVER
## 9 9 13 PASIR PANJANG 1 QTSZ13 N QUEENSTOWN
## 10 10 7 QUEENSWAY QTSZ07 N QUEENSTOWN
## PLN_AREA_C REGION_N REGION_C INC_CRC FMEL_UPD_D X_ADDR
## 1 MS CENTRAL REGION CR 5ED7EB253F99252E 2014-12-05 31595.84
## 2 OT CENTRAL REGION CR 8C7149B9EB32EEFC 2014-12-05 28679.06
## 3 SR CENTRAL REGION CR C35FEFF02B13E0E5 2014-12-05 29654.96
## 4 BM CENTRAL REGION CR 3775D82C5DDBEFBD 2014-12-05 26782.83
## 5 BM CENTRAL REGION CR 85D9ABEF0A40678F 2014-12-05 26201.96
## 6 BM CENTRAL REGION CR 9D286521EF5E3B59 2014-12-05 25358.82
## 7 BM CENTRAL REGION CR 7839A8577144EFE2 2014-12-05 27680.06
## 8 SR CENTRAL REGION CR 48661DC0FBA09F7A 2014-12-05 29253.21
## 9 QT CENTRAL REGION CR 1F721290C421BFAB 2014-12-05 22077.34
## 10 QT CENTRAL REGION CR 3580D2AFFBEE914C 2014-12-05 24168.31
## Y_ADDR SHAPE_Leng SHAPE_Area PA Child ECONOMY ACTIVE AGED
## 1 29220.19 5267.381 1630379.3 <NA> NA NA NA
## 2 29782.05 3506.107 559816.2 OUTRAM 260 4160 1280
## 3 29974.66 1740.926 160807.5 SINGAPORE RIVER 0 70 0
## 4 29933.77 3313.625 595428.9 BUKIT MERAH 430 6140 1690
## 5 30005.70 2825.594 387429.4 BUKIT MERAH 640 6580 1600
## 6 29991.38 4428.913 1030378.8 BUKIT MERAH 510 7600 1880
## 7 30230.86 3275.312 551732.0 BUKIT MERAH 670 8270 2230
## 8 30222.86 2208.619 290184.7 SINGAPORE RIVER 0 40 0
## 9 29893.78 6571.323 1084792.3 QUEENSTOWN 190 2270 330
## 10 30104.18 3454.239 631644.3 QUEENSTOWN 0 170 0
## TOTAL DEPENDENCY geometry
## 1 NA NA MULTIPOLYGON (((31495.56 30...
## 2 7510 0.3701923 MULTIPOLYGON (((29092.28 30...
## 3 80 0.0000000 MULTIPOLYGON (((29932.33 29...
## 4 11430 0.3452769 MULTIPOLYGON (((27131.28 30...
## 5 11760 0.3404255 MULTIPOLYGON (((26451.03 30...
## 6 14310 0.3144737 MULTIPOLYGON (((25899.7 297...
## 7 15180 0.3506651 MULTIPOLYGON (((27746.95 30...
## 8 50 0.0000000 MULTIPOLYGON (((29351.26 29...
## 9 4000 0.2290749 MULTIPOLYGON (((20996.49 30...
## 10 270 0.0000000 MULTIPOLYGON (((24472.11 29...Analysing Spatial Point Process
let’s understand the dataset ### 2017
summary(mpszpop2017$Child)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 120.0 450.0 806.7 1140.0 5710.02019
summary(mpszpop2019$Child)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 130.0 415.0 793.1 1065.0 5750.0####Linear Model
lm_2019 <- lm(Child ~TOTAL, data = mpsz2019)
lm_2019##
## Call:
## lm(formula = Child ~ TOTAL, data = mpsz2019)
##
## Coefficients:
## (Intercept) TOTAL
## -40.33533 0.05098###2017
ggplot(mpsz2019, aes(x=TOTAL, y=Child))+
geom_point()
####Linear Model
lm_2017 <- lm(Child ~TOTAL, data = mpsz2017)
lm_2017##
## Call:
## lm(formula = Child ~ TOTAL, data = mpsz2017)
##
## Coefficients:
## (Intercept) TOTAL
## -12.76578 0.05069ggplot(mpsz2017, aes(x=TOTAL, y=Child))+
geom_point()
### 2017 Childcare planning by subzone
Clean the childcare data set
Import childcare center datasets
childcare_2017 <- st_read(dsn = "data/2017/geospatial", layer="CHILDCARE")## Reading layer `CHILDCARE' from data source `C:\Users\sunho\Documents\SMU\3.2\gis\Take home\Take_Home\data\2017\geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 1312 features and 18 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 11203.01 ymin: 25667.6 xmax: 45404.24 ymax: 49300.88
## projected CRS: SVY21childcare_2019 <- st_read(dsn ="data/2020/Geospatial", layer="Childcare")## Reading layer `Childcare' from data source `C:\Users\sunho\Documents\SMU\3.2\gis\Take home\Take_Home\data\2020\Geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 1545 features and 25 fields
## geometry type: POINT
## dimension: XYZ
## bbox: xmin: 11203.01 ymin: 25667.6 xmax: 45404.24 ymax: 49300.88
## z_range: zmin: 0 zmax: 0
## projected CRS: SVY21 / Singapore TMCheck for CRS
let’s check if the CRS is correct format
st_crs(childcare_2017)## Coordinate Reference System:
## User input: SVY21
## wkt:
## PROJCRS["SVY21",
## BASEGEOGCRS["SVY21[WGS84]",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]st_crs(childcare_2019)## Coordinate Reference System:
## User input: SVY21 / Singapore TM
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]As the EPSG is 9001, convert it into 3414
childcare_2017 <- st_set_crs(childcare_2017, 3414)
childcare_2019 <- st_set_crs(childcare_2019, 3414)Check for CRS
Check once again to make sure CRS is correct format
st_crs(childcare_2017)## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]st_crs(childcare_2019)## Coordinate Reference System:
## User input: EPSG:3414
## wkt:
## PROJCRS["SVY21 / Singapore TM",
## BASEGEOGCRS["SVY21",
## DATUM["SVY21",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4757]],
## CONVERSION["Singapore Transverse Mercator",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",1.36666666666667,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",103.833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",28001.642,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",38744.572,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (N)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (E)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["Singapore"],
## BBOX[1.13,103.59,1.47,104.07]],
## ID["EPSG",3414]]filter the data
Filter the data to what is needed for analysis
childcare_2017 <- childcare_2017 %>% dplyr::select(ADDRESSPOS, ADDRESSSTR, NAME, geometry)
childcare_2019 <- childcare_2019 %>% dplyr::select(ADDRESSPOS, ADDRESSSTR, Name, geometry)check if there is NA data
Make sure that there is no NA data.
childcare_2017[rowSums(is.na(childcare_2017))!=0,]## Simple feature collection with 0 features and 3 fields
## bbox: xmin: NA ymin: NA xmax: NA ymax: NA
## projected CRS: SVY21 / Singapore TM
## [1] ADDRESSPOS ADDRESSSTR NAME geometry
## <0 rows> (or 0-length row.names)childcare_2019[rowSums(is.na(childcare_2019))!=0,]## Simple feature collection with 0 features and 3 fields
## bbox: xmin: NA ymin: NA xmax: NA ymax: NA
## projected CRS: SVY21 / Singapore TM
## [1] ADDRESSPOS ADDRESSSTR Name geometry
## <0 rows> (or 0-length row.names)Check for dublicated from the dataset
See if there is duplicated data as this can mislead the result
str(childcare_2017[duplicated(childcare_2017[,c('NAME')]), ]) ## Classes 'sf' and 'data.frame': 0 obs. of 4 variables:
## $ ADDRESSPOS: chr
## $ ADDRESSSTR: chr
## $ NAME : chr
## $ geometry :sfc_GEOMETRY of length 0 - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA
## ..- attr(*, "names")= chr [1:3] "ADDRESSPOS" "ADDRESSSTR" "NAME"str(childcare_2019[duplicated(childcare_2019[,c('Name')]), ]) ## Classes 'sf' and 'data.frame': 0 obs. of 4 variables:
## $ ADDRESSPOS: chr
## $ ADDRESSSTR: chr
## $ Name : chr
## $ geometry :sfc_GEOMETRY of length 0 - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA
## ..- attr(*, "names")= chr [1:3] "ADDRESSPOS" "ADDRESSSTR" "Name"Plot the childcare to understand which regions have childcare centres
2017
plot(childcare_2017)
#### 2019
plot(childcare_2019)
#Demand and Supply there was 1312 childcare centers in 2017
summary(childcare_2017)## ADDRESSPOS ADDRESSSTR NAME geometry
## Length:1312 Length:1312 Length:1312 POINT :1312
## Class :character Class :character Class :character epsg:3414 : 0
## Mode :character Mode :character Mode :character +proj=tmer...: 0There is 1545 childcare centers in 2019
the childcare centers have increased by 233 this means that there are higher demand in 2019 compare to 2017.
summary(childcare_2019)## ADDRESSPOS ADDRESSSTR Name geometry
## Length:1545 Length:1545 Length:1545 POINT Z :1545
## Class :character Class :character Class :character epsg:3414 : 0
## Mode :character Mode :character Mode :character +proj=tmer...: 0plotting choropleth maps 2017
tm_shape(mpsz2017)+
tm_fill("Child",
n = 8,
style = "jenks") +
tm_borders(alpha = 0.5)
#### plotting choropleth maps 2019
tm_shape(mpsz2019)+
tm_fill("Child",
n = 8,
style = "jenks") +
tm_borders(alpha = 0.5)
If I compare these two maps of child population distribution, I can understand that lesser were born in East, North-east and North. These means That childcare should be opened more in those areas
Section B
Convert sf dataset into ppp format
In order to do exploratory spatial data analysis, it is important to convert the sf dataset into ppp format.
firstly, convert the sf dataframe into SpatialPoints dataframe Secondly, SpatialPoints to generic sp format Lastly, generic sp to spatstat’s ppp format
codes are shown below
childcare_2017_sf <- as_Spatial(childcare_2017)
childcare_2017_sp <- as(childcare_2017_sf, "SpatialPoints")
childcare_ppp2017 <- as(childcare_2017_sp, "ppp")
childcare_2019_sf <- as_Spatial(childcare_2019)
childcare_2019_sp <- as(childcare_2019_sf, "SpatialPoints")
childcare_ppp2019 <- as(childcare_2019_sp, "ppp")Clean the duplicated points
2017
check if the childcare_ppp dataset have duplicated
any(duplicated(childcare_ppp2017))## [1] TRUE2019
any(duplicated(childcare_ppp2019))## [1] TRUEAs both 2017 and 2019 have the duplicated point. Let’s find out how many duplicated points are in 2017 and 2019 ### number of duplicated in 2017
sum(multiplicity(childcare_ppp2017) > 1)## [1] 852017 ppp dataset have 85 duplicated data
number of duplicated in 2017
sum(multiplicity(childcare_ppp2019) > 1)## [1] 1282019 ppp dataset have 128 duplicated data
Remove the duplicated points
Jittering approach
I will be using jilttering method to overcome the overlap issue
childcare_ppp_jit_2017 <- rjitter(childcare_ppp2017, retry=TRUE, nsim=1, drop=TRUE)
childcare_ppp_jit_2019 <- rjitter(childcare_ppp2019, retry=TRUE, nsim=1, drop=TRUE)Check again if there is duplicated data
sum(multiplicity(childcare_ppp_jit_2017) > 1)## [1] 0number of duplicated in 2017
sum(multiplicity(childcare_ppp_jit_2019) > 1)## [1] 0plot the map of no duplicated points
par(mfrow=c(1,2))
plot(childcare_ppp_jit_2017, main = "2017")
plot(childcare_ppp_jit_2019, main = "2019")
It is hard to identify the regions without the Singapore layer
Combine singapore map and Childcare
Import that sg coastal line Spatstat requires the analytical data in ppp object form. Hence,convert a SpatialDataFrame into ppp object. We need to convert the SpatialDataFrame into Spatial object first.
sg <- readOGR(dsn = "data/2017/geospatial", layer="CostalOutline")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\sunho\Documents\SMU\3.2\gis\Take home\Take_Home\data\2017\geospatial", layer: "CostalOutline"
## with 60 features
## It has 4 fieldssg_sp <- as(sg, "SpatialPolygons")Convert generic sp format to owin objection to allow singapore layer and childcare ppp to be combined
sg_owin <- as(sg_sp,"owin")Combind Childcare points with the coataloutline
childcareSG_ppp_2017 = childcare_ppp_jit_2017[sg_owin]
childcareSG_ppp_2019 = childcare_ppp_jit_2019[sg_owin]Plot the map
Now we can identify where the childcare centeres are located in Singapre
par(mfrow=c(1,2))
plot(childcareSG_ppp_2017, main = "2017")
plot(childcareSG_ppp_2019, main = "2019")
Chi-Square test
2017
qt_2017 <- quadrat.test(childcareSG_ppp_2017,
nx = 20, ny = 15)
qt_2017##
## Chi-squared test of CSR using quadrat counts
##
## data: childcareSG_ppp_2017
## X2 = 2027.7, df = 184, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)2019
qt_2019 <- quadrat.test(childcareSG_ppp_2019,
nx = 20, ny = 15)
qt_2019##
## Chi-squared test of CSR using quadrat counts
##
## data: childcareSG_ppp_2019
## X2 = 2433, df = 184, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 185 tiles (irregular windows)par(mfrow=c(1,2))
plot(childcareSG_ppp_2017)
plot(qt_2017, add = TRUE, cex =.1)
plot(childcareSG_ppp_2019)
plot(qt_2019, add = TRUE, cex =.1)
# Hyphothesis test
Using G-Function
This allow me to measures the distribution of the distances from an arbitary events to it’s nearest event.
par(mfrow=c(1,2))
G_2017 = Gest(childcareSG_ppp_2017, correction = "border")
plot(G_2017)
G_2019 = Gest(childcareSG_ppp_2019, correction = "border")
plot(G_2019)
to understand the spatial patterns, I will conduct hyphothesis test.
2017
H0 = the distribution of childcare service is randomly distributed in 2017 H1 = the distribution of childcare service is randomly distributed in 2017
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
G_2017.csr <- envelope(childcareSG_ppp_2017, Gest, nsim = 999)## Generating 999 simulations of CSR ...
## 1, 2, 3, ......10 [etd 3:43] .........20 [etd 3:42] .........
## 30 [etd 3:39] .........40 [etd 3:36] .........50 [etd 3:34] ........
## .60 [etd 3:31] .........70 [etd 3:30] .........80 [etd 3:28] .......
## ..90 [etd 3:26] .........100 [etd 3:24] .........110 [etd 3:22] ......
## ...120 [etd 3:19] .........130 [etd 3:17] .........140 [etd 3:15] .....
## ....150 [etd 3:14] .........160 [etd 3:11] .........170 [etd 3:09] ....
## .....180 [etd 3:06] .........190 [etd 3:03] .........200 [etd 3:01] ...
## ......210 [etd 2:59] .........220 [etd 2:56] .........230 [etd 2:54] ..
## .......240 [etd 2:51] .........250 [etd 2:49] .........260 [etd 2:47] .
## ........270 [etd 2:44] .........280 [etd 2:42] .........290
## [etd 2:39] .........300 [etd 2:37] .........310 [etd 2:35] .........
## 320 [etd 2:32] .........330 [etd 2:30] .........340 [etd 2:28] ........
## .350 [etd 2:26] .........360 [etd 2:23] .........370 [etd 2:21] .......
## ..380 [etd 2:19] .........390 [etd 2:17] .........400 [etd 2:14] ......
## ...410 [etd 2:12] .........420 [etd 2:10] .........430 [etd 2:07] .....
## ....440 [etd 2:05] .........450 [etd 2:03] .........460 [etd 2:01] ....
## .....470 [etd 1:58] .........480 [etd 1:56] .........490 [etd 1:54] ...
## ......500 [etd 1:52] .........510 [etd 1:49] .........520 [etd 1:47] ..
## .......530 [etd 1:45] .........540 [etd 1:42] .........550 [etd 1:40] .
## ........560 [etd 1:38] .........570 [etd 1:36] .........580
## [etd 1:34] .........590 [etd 1:31] .........600 [etd 1:29] .........
## 610 [etd 1:27] .........620 [etd 1:25] .........630 [etd 1:22] ........
## .640 [etd 1:20] .........650 [etd 1:18] .........660 [etd 1:16] .......
## ..670 [etd 1:13] .........680 [etd 1:11] .........690 [etd 1:09] ......
## ...700 [etd 1:07] .........710 [etd 1:05] .........720 [etd 1:02] .....
## ....730 [etd 1:00] .........740 [etd 58 sec] .........750 [etd 56 sec] ....
## .....760 [etd 53 sec] .........770 [etd 51 sec] .........780 [etd 49 sec] ...
## ......790 [etd 47 sec] .........800 [etd 44 sec] .........810 [etd 42 sec] ..
## .......820 [etd 40 sec] .........830 [etd 38 sec] .........840 [etd 36 sec] .
## ........850 [etd 33 sec] .........860 [etd 31 sec] .........870
## [etd 29 sec] .........880 [etd 27 sec] .........890 [etd 25 sec] .........
## 900 [etd 22 sec] .........910 [etd 20 sec] .........920 [etd 18 sec] ........
## .930 [etd 16 sec] .........940 [etd 13 sec] .........950 [etd 11 sec] .......
## ..960 [etd 9 sec] .........970 [etd 7 sec] .........980 [etd 4 sec] ......
## ...990 [etd 2 sec] ........ 999.
##
## Done.plot(G_2017.csr)
The point pattern of Government-linked childcare centres in 2017 largely lies within the envelope, which means it is not statistically significant and we fail to reject the Ho.
2019
H0 = the distribution of childcare service is randomly distributed in 2019 H1 = the distribution of childcare service is randomly distributed in 2019
The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.
G_2019.csr <- envelope(childcareSG_ppp_2019, Gest, nsim = 999)## Generating 999 simulations of CSR ...
## 1, 2, 3, ......10 [etd 4:14] .........20 [etd 4:10] .........
## 30 [etd 4:08] .........40 [etd 4:05] .........50 [etd 4:02] ........
## .60 [etd 4:00] .........70 [etd 3:57] .........80 [etd 3:55] .......
## ..90 [etd 3:52] .........100 [etd 3:49] .........110 [etd 3:47] ......
## ...120 [etd 3:44] .........130 [etd 3:42] .........140 [etd 3:41] .....
## ....150 [etd 3:38] .........160 [etd 3:35] .........170 [etd 3:32] ....
## .....180 [etd 3:30] .........190 [etd 3:27] .........200 [etd 3:24] ...
## ......210 [etd 3:22] .........220 [etd 3:19] .........230 [etd 3:17] ..
## .......240 [etd 3:14] .........250 [etd 3:11] .........260 [etd 3:09] .
## ........270 [etd 3:06] .........280 [etd 3:04] .........290
## [etd 3:01] .........300 [etd 2:59] .........310 [etd 2:57] .........
## 320 [etd 2:54] .........330 [etd 2:52] .........340 [etd 2:49] ........
## .350 [etd 2:46] .........360 [etd 2:44] .........370 [etd 2:41] .......
## ..380 [etd 2:39] .........390 [etd 2:36] .........400 [etd 2:33] ......
## ...410 [etd 2:31] .........420 [etd 2:28] .........430 [etd 2:26] .....
## ....440 [etd 2:23] .........450 [etd 2:21] .........460 [etd 2:18] ....
## .....470 [etd 2:15] .........480 [etd 2:13] .........490 [etd 2:11] ...
## ......500 [etd 2:08] .........510 [etd 2:05] .........520 [etd 2:03] ..
## .......530 [etd 2:00] .........540 [etd 1:58] .........550 [etd 1:55] .
## ........560 [etd 1:52] .........570 [etd 1:50] .........580
## [etd 1:47] .........590 [etd 1:45] .........600 [etd 1:42] .........
## 610 [etd 1:40] .........620 [etd 1:37] .........630 [etd 1:34] ........
## .640 [etd 1:32] .........650 [etd 1:29] .........660 [etd 1:27] .......
## ..670 [etd 1:24] .........680 [etd 1:22] .........690 [etd 1:19] ......
## ...700 [etd 1:17] .........710 [etd 1:14] .........720 [etd 1:11] .....
## ....730 [etd 1:09] .........740 [etd 1:06] .........750 [etd 1:04] ....
## .....760 [etd 1:01] .........770 [etd 59 sec] .........780 [etd 56 sec] ...
## ......790 [etd 54 sec] .........800 [etd 51 sec] .........810 [etd 48 sec] ..
## .......820 [etd 46 sec] .........830 [etd 43 sec] .........840 [etd 41 sec] .
## ........850 [etd 38 sec] .........860 [etd 36 sec] .........870
## [etd 33 sec] .........880 [etd 31 sec] .........890 [etd 28 sec] .........
## 900 [etd 25 sec] .........910 [etd 23 sec] .........920 [etd 20 sec] ........
## .930 [etd 18 sec] .........940 [etd 15 sec] .........950 [etd 13 sec] .......
## ..960 [etd 10 sec] .........970 [etd 7 sec] .........980 [etd 5 sec] ......
## ...990 [etd 2 sec] ........ 999.
##
## Done.plot(G_2019.csr)
The point pattern of Government-linked childcare centres in 2019 largely lies within the envelope, which means it is not statistically significant and we fail to reject the Ho.
Kernel Density Estimation
Computing kernel density estimation
I am doing to use automatic bandwidth selection method
bw_2017 <- bw.diggle(childcareSG_ppp_2017)
bw_2019 <- bw.diggle(childcareSG_ppp_2019)kde_childcareSG_bw_2017 <- density(childcareSG_ppp_2017, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
kde_childcareSG_bw_2019 <- density(childcareSG_ppp_2019, sigma=bw.diggle, edge=TRUE, kernel="gaussian") Plot to display the kernel density
More childcare centers are opened since 2017 especially at East, North-east and North.
par(mfrow=c(1,2))
plot(kde_childcareSG_bw_2017)
plot(kde_childcareSG_bw_2019)
However, range is too small to comprehended hence let’s rerange the scale.
childcareSG_ppp_2017.km <- rescale(childcareSG_ppp_2017, 1000, "km")
childcareSG_ppp_2019.km <- rescale(childcareSG_ppp_2019, 1000, "km")par(mfrow=c(1,2))
kde_childcareSG_2017.bw_2017 <- density(childcareSG_ppp_2017.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_2017.bw_2017)
kde_childcareSG_2019.bw_2019 <- density(childcareSG_ppp_2019.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_2019.bw_2019)
Work on the automatic badwidth method
I am going to use bw.ppl() as it produce better appropriate values when the pattern consists predominatly of tight cluters
2017
bw.ppl(childcareSG_ppp_2017.km)## sigma
## 0.40021382019
bw.ppl(childcareSG_ppp_2019.km)## sigma
## 0.4136135Comput kernel density estimation
par(mfrow=c(1,2))
kde_childcareSG_600_2017 <- density(childcareSG_ppp_2017.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600_2017)
kde_childcareSG_600_2019 <- density(childcareSG_ppp_2019.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600_2019)
Convert KDE output into grid object
2017
gridded_kde_childcareSG_bw_2017 <- as.SpatialGridDataFrame.im(kde_childcareSG_2017.bw_2017)
spplot(gridded_kde_childcareSG_bw_2017)
2019
gridded_kde_childcareSG_bw_2019 <- as.SpatialGridDataFrame.im(kde_childcareSG_2019.bw_2019)
spplot(gridded_kde_childcareSG_bw_2019)
## Conver gridded output into raster
kde_childcareSG_bw_raster_2017 <- raster(gridded_kde_childcareSG_bw_2017)
kde_childcareSG_bw_raster_2019 <- raster(gridded_kde_childcareSG_bw_2019)let’s understand the properties of the RasterLayer 2017
kde_childcareSG_bw_raster_2017## class : RasterLayer
## dimensions : 128, 128, 16384 (nrow, ncol, ncell)
## resolution : 0.4170614, 0.2647348 (x, y)
## extent : 2.663926, 56.04779, 16.35798, 50.24403 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : v
## values : -4.388506e-15, 29.18322 (min, max)Crs is showing NA let’s add CRS information
projection(kde_childcareSG_bw_raster_2017) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster_2017## class : RasterLayer
## dimensions : 128, 128, 16384 (nrow, ncol, ncell)
## resolution : 0.4170614, 0.2647348 (x, y)
## extent : 2.663926, 56.04779, 16.35798, 50.24403 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
## source : memory
## names : v
## values : -4.388506e-15, 29.18322 (min, max)2019
kde_childcareSG_bw_raster_2019## class : RasterLayer
## dimensions : 128, 128, 16384 (nrow, ncol, ncell)
## resolution : 0.4170614, 0.2647348 (x, y)
## extent : 2.663926, 56.04779, 16.35798, 50.24403 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : v
## values : -7.51131e-15, 26.62348 (min, max)Crs is showing NA let’s add CRS information
projection(kde_childcareSG_bw_raster_2019) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster_2019## class : RasterLayer
## dimensions : 128, 128, 16384 (nrow, ncol, ncell)
## resolution : 0.4170614, 0.2647348 (x, y)
## extent : 2.663926, 56.04779, 16.35798, 50.24403 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
## source : memory
## names : v
## values : -7.51131e-15, 26.62348 (min, max)Visualise the tmap
2017
tm_shape(kde_childcareSG_bw_raster_2017) +
tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE, main.title = 2017)
tm_shape(kde_childcareSG_bw_raster_2019) +
tm_raster("v") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE, main.title = 2019)
Compare the advantage of kernel density maps over point maps
Advantage of Kernel Density Map
Firstly, Kernel density map allow me to understand the planning area with different intensity. compare to spatial map, if I want to identify which areas are more intense, I need to count the dots and it is hard to identify the exact locations unlike kernal map allow us to see the location clearly.
Furthermore, for point map, if the dots are lesser, it is easy to identify the location but once they have a lot of dots, it is difficult to understand the data.