Background

With the inception of The Smart Nation vision in 2017, Open Data is seen as a necessary component of this initiative, especially in the promotion of public-private collaborations (co-innovation).

The specific objectives of the study are as follow:

  1. To gain understanding on the supply and demand of childcare and kindergarten services at planning subzone level.

  2. To gain understanding on the geographic distribution of childcare and kindergarten in Singapore.

Packages Installation

packages = c('rgdal', 'maptools', 'raster','spatstat', 'tmap','sf', 'dplyr', 'tidyverse','units','ggplot2','OpenStreetMap')
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.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/sihua/OneDrive/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/sihua/OneDrive/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Loading required package: maptools
## Checking rgeos availability: TRUE
## Loading required package: raster
## Loading required package: spatstat
## Loading required package: spatstat.data
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## Loading required package: rpart
## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.64-1 is out of date by more than 4 months; we recommend upgrading to the latest version.
## 
## Attaching package: 'spatstat'
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
## Loading required package: tmap
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: tidyverse
## Registered S3 method overwritten by 'cli':
##   method     from    
##   print.boxx spatstat
## -- Attaching packages ------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v stringr 1.4.0
## v tidyr   1.1.1     v forcats 0.5.0
## v readr   1.3.1
## -- Conflicts ---------------------------------------- tidyverse_conflicts() --
## x dplyr::collapse() masks nlme::collapse()
## x tidyr::extract()  masks raster::extract()
## x dplyr::filter()   masks stats::filter()
## x dplyr::lag()      masks stats::lag()
## x dplyr::select()   masks raster::select()
## Loading required package: units
## udunits system database from C:/Users/sihua/OneDrive/Documents/R/win-library/4.0/units/share/udunits
## Loading required package: OpenStreetMap

Geospatial Data Wrangling

Data Sources

This study will be using the following datasets:

Data Import

URA Master Plan 2014 Planning Subzone Dataset:

mpsz <- st_read(dsn =  "data/geospatial",layer="MP14_SUBZONE_WEB_PL")
## Reading layer `MP14_SUBZONE_WEB_PL' from data source `C:\GIS415\TakeHome\IS415_Take-home_Ex01_PengSzuHua\data\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:  SVY21
mpsz1<- readOGR(dsn =  "data/geospatial",layer="MP14_SUBZONE_WEB_PL")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\GIS415\TakeHome\IS415_Take-home_Ex01_PengSzuHua\data\geospatial", layer: "MP14_SUBZONE_WEB_PL"
## with 323 features
## It has 15 fields

ChildCare Dataset:

sf_childcare = st_read("data/geospatial/child-care-services-kml.kml")
## Reading layer `CHILDCARE' from data source `C:\GIS415\TakeHome\IS415_Take-home_Ex01_PengSzuHua\data\geospatial\child-care-services-kml.kml' using driver `KML'
## Simple feature collection with 1545 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

Kindergarten Dataset:

sf_kindergarten = st_read("data/geospatial/kindergartens-kml.kml")
## Reading layer `KINDERGARTENS' from data source `C:\GIS415\TakeHome\IS415_Take-home_Ex01_PengSzuHua\data\geospatial\kindergartens-kml.kml' using driver `KML'
## Simple feature collection with 448 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6887 ymin: 1.247759 xmax: 103.9717 ymax: 1.455452
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

Singapore Residents by Planning AreaSubzone, Age Group, Sex and Type of Dwelling Dataset:

popresplan <- read_csv("data/aspatial/respopagesextod2011to2019.csv")
## Parsed with column specification:
## cols(
##   PA = col_character(),
##   SZ = col_character(),
##   AG = col_character(),
##   Sex = col_character(),
##   TOD = col_character(),
##   Pop = col_double(),
##   Time = col_double()
## )

Working With Projection

We can see from the above result, all the datatset projection is WGS84. Hence,this section is to ensure that all the datasets are in the same projection system and EPSG. In this analysis, to aligned with the SG projection system we will be using SVY21:3414.

Subzone Dataset

Check Projection-Subzone DataSet

st_crs(mpsz1)
## 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]]]]

Projection Transformation-subzone

Using st_transform () function to transform the projection system from WGS84 to 3414.

Changing Projection into 3414 to mpsz simple feature dataframe and spaital point.

sp_mpsz3414 <- spTransform(mpsz1,CRS("+init=epsg:3414"))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum SVY21 in CRS definition
sf_mpsz3414 <- st_transform(mpsz,3414)

Ensure the subzone are transformed into EPSG:3414.

st_crs(sf_mpsz3414)
## 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]]

ChildCare Dataset

Check Projection-Childcare DataSet

st_crs(sf_childcare)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Projection Transformation - ChildCare

To transform childcare data point into spatial point, as_spatial() function is used.

sf_childcare3414 <- st_transform(sf_childcare, 3414)
sp_child3414<- as_Spatial(sf_childcare3414)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

Check Transformed Projection- ChildCare

st_crs(sf_childcare3414)
## 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]]

Kindergarten Dataset

Check Projection- Kindergarten DataSet

st_crs(sf_kindergarten)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Projection Transformation - Kindergarten

sf_kindergarten3414 <- st_transform(sf_kindergarten, 3414)

sp_kindergarten3414<- as_Spatial(sf_kindergarten3414)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

Check Transformed Projection- Kindergarten

st_crs(sf_kindergarten3414)
## 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]]

Data Preparation

The data preparation is preformed to:

  1. Extracting 2018 records only.

    • By using year 2018 records, the age group 0 to 4 will be 2- 6 years old in 2020. which also fulfill the childcare and kindergarten age group.

  2. Group by Planning area and subzone.

  3. Sum up the total population by subzone.

  4. To derive the following data:

    • Planning Area
    • Subzone
    • Age Group: 0 - 4
    • Age Group: 5 - 9
    • Update the cell values of Planning area and Subzone to upper case.
popresplan2018 <- popresplan %>%
  filter(Time == 2018)%>%
  group_by(PA,SZ,AG)%>%
  summarize(`POP`=sum(`Pop`))%>%
  spread(AG,POP)%>%
  select(`PA`, `SZ`,`0_to_4`,`5_to_9`)%>%
  mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper)) %>%

  print(popresplan2018)
## `summarise()` regrouping output by 'PA', 'SZ' (override with `.groups` argument)
## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## # A tibble: 323 x 4
## # Groups:   PA, SZ [323]
##    PA         SZ                     `0_to_4` `5_to_9`
##    <chr>      <chr>                     <dbl>    <dbl>
##  1 ANG MO KIO ANG MO KIO TOWN CENTRE      180      290
##  2 ANG MO KIO CHENG SAN                  1040     1080
##  3 ANG MO KIO CHONG BOON                  890      890
##  4 ANG MO KIO KEBUN BAHRU                 710      870
##  5 ANG MO KIO SEMBAWANG HILLS             200      310
##  6 ANG MO KIO SHANGRI-LA                  560      640
##  7 ANG MO KIO TAGORE                      250      330
##  8 ANG MO KIO TOWNSVILLE                  820      950
##  9 ANG MO KIO YIO CHU KANG                  0        0
## 10 ANG MO KIO YIO CHU KANG EAST           160      150
## # ... with 313 more rows

From the above result, we can see that there are total of 323 subzone in Singapore.

Joining the attribute data and geospatial data

In order to clip the URA subzone map with the popresplan2019 data, left_join() function of dplyr is used to join the geographical data and attribute table using planning subzone.

mpsz_popresplan2018 <- left_join(sf_mpsz3414, popresplan2018, by = c("SUBZONE_N" = "SZ"))

print(mpsz_popresplan2018)
## Simple feature collection with 323 features and 18 fields
## geometry type:  MULTIPOLYGON
## 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 0_to_4 5_to_9
## 1  29220.19   5267.381  1630379.3    MARINA SOUTH      0      0
## 2  29782.05   3506.107   559816.2          OUTRAM    230    240
## 3  29974.66   1740.926   160807.5 SINGAPORE RIVER      0      0
## 4  29933.77   3313.625   595428.9     BUKIT MERAH    370    460
## 5  30005.70   2825.594   387429.4     BUKIT MERAH    590    790
## 6  29991.38   4428.913  1030378.8     BUKIT MERAH    460    490
## 7  30230.86   3275.312   551732.0     BUKIT MERAH    580    640
## 8  30222.86   2208.619   290184.7 SINGAPORE RIVER     10      0
## 9  29893.78   6571.323  1084792.3      QUEENSTOWN    180    190
## 10 30104.18   3454.239   631644.3      QUEENSTOWN      0      0
##                          geometry
## 1  MULTIPOLYGON (((31495.56 30...
## 2  MULTIPOLYGON (((29092.28 30...
## 3  MULTIPOLYGON (((29932.33 29...
## 4  MULTIPOLYGON (((27131.28 30...
## 5  MULTIPOLYGON (((26451.03 30...
## 6  MULTIPOLYGON (((25899.7 297...
## 7  MULTIPOLYGON (((27746.95 30...
## 8  MULTIPOLYGON (((29351.26 29...
## 9  MULTIPOLYGON (((20996.49 30...
## 10 MULTIPOLYGON (((24472.11 29...

Supply and demand of childcare and kindergarten services by planning subzone

Exploratory Spatial Data Analysis

This Section will be perform EDA on the supply and demand of childcare and kindergarten services by planning subzone.

Point In Polygon By Services

To see the distribution of the childcare and kindergarten services, It is required to understand the number of services per planning subzone.Thus, st_intersection() function will be used to count the number of service point in polygon by mapping the respective services geometry point with the planning subzone layer.

ChildCare Services

mpsz_popresplan2018$`childcare Count`<- lengths(st_intersects(mpsz_popresplan2018, sf_childcare3414))

Kindergarten Services

mpsz_popresplan2018$`Kindergarten Count`<- lengths(st_intersects(mpsz_popresplan2018, sf_kindergarten3414))
print(mpsz_popresplan2018)
## Simple feature collection with 323 features and 20 fields
## geometry type:  MULTIPOLYGON
## 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 0_to_4 5_to_9
## 1  29220.19   5267.381  1630379.3    MARINA SOUTH      0      0
## 2  29782.05   3506.107   559816.2          OUTRAM    230    240
## 3  29974.66   1740.926   160807.5 SINGAPORE RIVER      0      0
## 4  29933.77   3313.625   595428.9     BUKIT MERAH    370    460
## 5  30005.70   2825.594   387429.4     BUKIT MERAH    590    790
## 6  29991.38   4428.913  1030378.8     BUKIT MERAH    460    490
## 7  30230.86   3275.312   551732.0     BUKIT MERAH    580    640
## 8  30222.86   2208.619   290184.7 SINGAPORE RIVER     10      0
## 9  29893.78   6571.323  1084792.3      QUEENSTOWN    180    190
## 10 30104.18   3454.239   631644.3      QUEENSTOWN      0      0
##                          geometry childcare Count Kindergarten Count
## 1  MULTIPOLYGON (((31495.56 30...               0                  0
## 2  MULTIPOLYGON (((29092.28 30...               4                  1
## 3  MULTIPOLYGON (((29932.33 29...               0                  0
## 4  MULTIPOLYGON (((27131.28 30...               3                  1
## 5  MULTIPOLYGON (((26451.03 30...               1                  1
## 6  MULTIPOLYGON (((25899.7 297...               8                  2
## 7  MULTIPOLYGON (((27746.95 30...               3                  1
## 8  MULTIPOLYGON (((29351.26 29...               2                  0
## 9  MULTIPOLYGON (((20996.49 30...               2                  5
## 10 MULTIPOLYGON (((24472.11 29...               1                  0

Calculate Land Area Per Subzone

mpsz_popresplan2018$LandArea <- mpsz_popresplan2018 %>%
  st_area()%>%
  set_units(km^2)

Calculate Density

mpsz_popresplan2018 <- mpsz_popresplan2018 %>%
mutate(
      `Childcare Density` = as.numeric(`childcare Count`/LandArea),
      `Kindergarten Density` = as.numeric(`Kindergarten Count`/LandArea)
      )

Replacing 0 values with NA

It is required to replace 0 values with NA as the default missing value will be shaded in gray color and we will be able to differentiate 0 values from the lower range values (example : values with 1).

mpsz_popresplan2018 <- mpsz_popresplan2018 %>%
   mutate(
          `childcare Count` = na_if(`childcare Count`,0),
          `Kindergarten Count` = na_if(`Kindergarten Count`,0),
          `Childcare Density` = na_if(`Childcare Density`,0),
          `Kindergarten Density` = na_if(`Kindergarten Density`,0),
          `0_to_4` = na_if(`0_to_4`,0),
          `5_to_9` = na_if(`5_to_9`,0),
          `childcareRatio` = as.numeric(`0_to_4`/`childcare Count`),
          `KindergartenRatio` = as.numeric(`0_to_4`/`Kindergarten Count`)

        )

Choropleth Map by Number of Services

In this section choropleth map will be used to display the distribution on the number of childcare and kindergarten services and also density of the services per subzone.

ChildCare Services

Distribution on Number of childcare services

Childcaredistribution<-tm_shape(mpsz_popresplan2018)+
  tm_fill("childcare Count",
          breaks = c(1, 11, 21, 31, 41,50))+
  tm_borders(alpha = 0.5)+
  tm_layout(main.title = "Distribution Map: # of childcare services",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE)

Density on ChildCare Services:

ChildCareDensitymap<-tm_shape(mpsz_popresplan2018)+
  tm_polygons("Childcare Density")+
  tm_layout(main.title = "Density Map: # of Childcare Services",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
           frame = TRUE)

Kindergarten Services

Distribution on Number of Kindergarten services

Kindergartendistribution<-tm_shape(mpsz_popresplan2018)+
  tm_fill("Kindergarten Count",
          breaks = c(1, 11, 21, 31, 41,50))+
  tm_borders(alpha = 0.5)+
  tm_layout(main.title = "Distribution Map: # of Kindergarten Services",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
           frame = TRUE)

Density on Kindergarten Services:

KindergartenDensitymap<-tm_shape(mpsz_popresplan2018)+
  tm_polygons("Kindergarten Density")+
  tm_layout(main.title = "Density Map: # of Kindergarten Services",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
           frame = TRUE)

Let’s look at the overview on the Density map by number of services. To have a better overview, i have bin the number of services on to the following range: “1-10”, “11-20”,“21-30”,“31-40” and “41-50” so that we will be able to make comparison on the same scale.

Primarily Service Distribution Overview:

  • Distribution of both services are similar, however child care services are more wide spread across Singapore compared to Kindergarten services.

  • Generally, the number of child care services has a higher density in different subzone compared to kindergarten services. As you can see from the density map, the higherest number of service in subzone for child care services are in the range of 41 to 50 however,the highest number of kindergarten services only in range 11 to 20.

tmap_arrange(Childcaredistribution, Kindergartendistribution, ChildCareDensitymap,KindergartenDensitymap)

Based on the current services distribution and the density of the services per subzone (based on number of service count) pattern , we can see that Child care services have more supply compared to Kindergarten services.

Choropleth Map by Population Age Group

To derive demand of the services, we can look at population. As the age group population will be able to determine the demand of the services.

Age 0 to 4

Distribution on age group 0 to 4:

Distribution0_to_4<-tm_shape(mpsz_popresplan2018)+
   tm_fill("0_to_4",
          breaks = c(1, 1001, 2001, 3001, 4001,5001,6000))+
  tm_borders(alpha = 0.5)+
  tm_layout(main.title = "Distibution Map:# of Age Group 0-4",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE)

Let’s look at the overview on the Density map by 2 population age group which are 1) 0 to 4 years old and 2)5 to 9 years old. To have a better overview, i have bin the number of population into the following range: “1-1000”, “1001-2000”,“2001-3000”,“3001-4000” and “4001-5000” and “5001-6000” so that we will be able to make comparison on the same scale.

Primarily Age Population Distribution Overview:

  • Based on the distribution map, both the distribution of age groups are having the same pattern.

  • The number of the population for these 2 age group are also sharing the same range. As you can see the for both highest number of the population age group are within 5001 to 6000.

tmap_arrange(Distribution0_to_4, asp=1, ncol=2)

Population Vs Number of services provided

Bar Plot On # of ChildCare Services by subzone

To plot a bar chart ggplot geom_bar() function is used.

bartchildcaresservices<-top_n(mpsz_popresplan2018, n=3,`childcare Count`) %>%
ggplot(., 
       aes(x=reorder(SUBZONE_N,-`childcare Count`), 
           y=`childcare Count`))+
geom_bar(stat='identity', 
         fill="steelblue")+
geom_text(aes(y=`childcare Count`, label=`childcare Count`), vjust=1.6, 
            color="white", size=3.5)+
xlab("Subzone")+
ylab("# of ChildCare Services")+
scale_fill_brewer(palette="Reds")+
theme_minimal()

Bar Plot On # of Kindergarten Services by subzone

bartkindergartenservices<-top_n(mpsz_popresplan2018, n= 3,`Kindergarten Count`) %>%
ggplot(., 
       aes(x=reorder(SUBZONE_N,-`Kindergarten Count`), 
           y=`Kindergarten Count`))+
geom_bar(stat='identity', 
         fill="steelblue")+
geom_text(aes(y=`Kindergarten Count`, label=`Kindergarten Count`), vjust=1.6, color="white", size=3.5)+
xlab("Subzone")+
ylab("# of Kindergarten Services")+
scale_fill_brewer(palette="Greens")

Bar Plot On # of Age 0 to 4 population by subzone

bar0to4<-top_n(mpsz_popresplan2018, n=3,`0_to_4`) %>%
ggplot(., 
       aes(x=reorder(SUBZONE_N,-`0_to_4`), 
           y=`0_to_4`))+
geom_bar(stat='identity', 
         fill="steelblue")+
geom_text(aes(y=`0_to_4`, label=`0_to_4`), vjust=1.6, 
            color="white", size=3.5)+
xlab("Subzone")+
ylab("Age 0 to 4  Poluation")

Analysis and Findings:

bar0to4

par(mfrow=c(1,2))
plot(bartchildcaresservices)

plot(bartkindergartenservices)

Based on the above Bar charts, it shows that the top 3 subzone that has the most number of Kindergarten services are

  • Tampines East -> 17 Kindergarten services
  • Katong -> 11 Kindergarten services
  • Aljunied , Frankel & Pasir Ris West -> 10 Kindergarten services

and the top 3 subzone that has the most childCare services are

  • Tampines East -> 42 childcare services
  • Woodlands East -> 42 childcare services
  • SengKang Town Centre -> 27 childcare services

and also the top 3 area of the highest population on age 2 - 6 years old are

  • Matilda ->5750
  • Waterway East -> 5730
  • Tampines East -> 5720

From the above result we can see that, Tampines East has the highest number of child care and kindergarten services is provided in the area and the population also in the top 3 list.

[Note: Population number does not have much difference in the number of population. Hence, those 3 area has high number of population age group from 2 - 6 years old.]

From the above result, it only shows that the supply and demand seems to fulfill in the Tampines East area only.However, is that true?

Supply Vs Demand on ChildCare and Kindergarten Services

Expected Ratio per service

Let’s look the ratio per service provider. To computed the ratio per service. capacity of the service provider is used. Based on ECDA statistic report as of April 2019, the total child care services capacity is 171,600 and total kindergarten capacity is 58,240.

The formula : $ Expected Ratio per Service = Total Service Capacity / Total Count of Service $

mpsz_popresplan2018ratio <- mpsz_popresplan2018 %>%
   mutate(
          `childcare Count` = replace_na(`childcare Count`,0),
          `Kindergarten Count` = replace_na(`Kindergarten Count`,0)
           )
Avgcapacitychildcare <- as.numeric(171600/sum(mpsz_popresplan2018ratio$`childcare Count`))
Avgcapacitykindergarten <- as.numeric(58240/sum(mpsz_popresplan2018ratio$`Kindergarten Count`))

Ratio per service by subzone

After computed the average capacity on the services, in order to determine if the service is in demand or over supply in that particular subzone the following formula will be used:

$ Ratio Per Service by Subzone = Age Group / Average Capacity$

mpsz_popresplan2018ratio <- mpsz_popresplan2018 %>%
   mutate(
          `childcareRatio` = as.numeric(`0_to_4`/Avgcapacitychildcare),
      `KindergartenRatio` = as.numeric(`0_to_4`/Avgcapacitykindergarten)
           )

Choropleth Map on Ratio Per Services

Supply and Demand on ChildCare Services by Subzone

ChildCareRatioPerServices<-tm_shape(mpsz_popresplan2018)+
   tm_fill("childcareRatio")+
  tm_borders(alpha = 0.5)+
  tm_layout(main.title = "ChildCare Ratio Per Service",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE)

Supply and Demand on Kindergarten Services by Subzone

KindergartenRatioPerServices<-tm_shape(mpsz_popresplan2018)+
   tm_fill("KindergartenRatio")+
  tm_borders(alpha = 0.5)+
  tm_layout(main.title = "Kindergarten Ratio Per Service",
            main.title.position = "center",
            main.title.size = 1.2,
            legend.height = 0.45, 
            legend.width = 0.35,
            frame = TRUE)

Supply & Demand Analysis

tmap_arrange(ChildCareRatioPerServices,KindergartenRatioPerServices, asp=1, ncol=2)

Conclusion :

  • Based on the population and distribution map, it shows that both services are well covered to the area with the necessary area.

  • Based on the Ratio chart it shows that the darker the color indicate that the subzone are in demand of the service and the lighter the color indicate that the subzone are over supply of the service.

  • Based on the Choropleth map, we can see that in overall childcare services are in higher demand compared to kindergarten services.

Spatial Point Pattern Analysis

In this section, the analysis is fosucing on spatial point pattern analysis with various satistical technics.

Exploratory Spatial Data Analysis

In this section, it will display the location of childcare and kindergarten services in separate point maps. Before plotting the data point, it is required to convert sf into sp and convert it to ppp object again.

sp_child<- as_Spatial(sf_childcare3414$geometry)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved
child_PPP<-as(spTransform(sp_child,CRS("+init=epsg:3414")), "ppp")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum SVY21 in CRS definition
sp_kindergarten<- as_Spatial(sf_kindergarten$geometry)
Kindergarten_PPP<-as(spTransform(sp_kindergarten,CRS("+init=epsg:3414")), "ppp")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum SVY21 in CRS definition
par(mfrow=c(1,2))
plot(child_PPP)
plot(Kindergarten_PPP)

## Handling duplicated points

This section will be checking if there is any duplicate point in the data.

any(duplicated(child_PPP))
## [1] TRUE
any(duplicated(Kindergarten_PPP))
## [1] TRUE

There are lots of way to handle the duplicate data,the easiest way is to delete the duplicates. However, at the same time the analysis might lost lots of useful point events.

Jittering is another way to handle duplicate, the function will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

childcare_ppp_jit <- rjitter(child_PPP, retry=TRUE, nsim=1, drop=TRUE)
Kindergarten_PPP_jit <- rjitter(Kindergarten_PPP, retry=TRUE, nsim=1, drop=TRUE)

Create Owins

When analysing spatial point patterns, Owin is created to represent polygonal region. It is a good practice to confine the analysis with a geographical area like Singapore boundary.

sg_owin <- as(sp_mpsz3414, "owin")
plot(sg_owin)

Combining service point with SG Polygon

This section will be combining the childcare service point and kindergarten service point with Singpore Polygon.

childcareSG_ppp = childcare_ppp_jit[sg_owin]
KindergartenSG_ppp = Kindergarten_PPP_jit[sg_owin]

par(mfrow=c(1,2))
plot(childcareSG_ppp)
plot(KindergartenSG_ppp)

From the map above, it shows that both child care services and kindergarten services are clustered at the same region as well as the population region that show in the distribution map in the earlier section. Thus, this shows that, the services provider are most likely located around the residential area.

Stasticstical Analysis

To confirm on the observation above, Conditional Monte Carlo test is used.

Monte Carlo test

This test is used to perform the quadrat analysis and compute variance-Mean Ratio (VMR) to indicate the type of the distribution.

Child Care Services

quadrat.test(childcareSG_ppp, 
             nx = 20, ny = 15,
             method="M",
             nsim=199)
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  childcareSG_ppp
## X2 = 2623.5, p-value = 0.01
## alternative hypothesis: two.sided
## 
## Quadrats: 193 tiles (irregular windows)

Based on the return X2= 2616.6 values we can concluded that childcare services is a clustered distribution as the variance is relatively large.

Kindergarten Services

quadrat.test(KindergartenSG_ppp, 
             nx = 20, ny = 15,
             method="M",
             nsim=999)
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  KindergartenSG_ppp
## X2 = 836.84, p-value = 0.004
## alternative hypothesis: two.sided
## 
## Quadrats: 193 tiles (irregular windows)

After performing conditional Monte Carlo test of CSR through 1000 simulations, the results remain consistent at rejecting the null hypothesis that the distribution of points are random.

Based on the return X2 values:895.58, we can concluded that kindergarten services is a strong clustered distribution as the variance is relatively large. However, kindergarten services is not as clustered as childcare services.

Nearest Neighbour Analysis

In this section, we will perform the Clark-Evans test of aggregation for a spatial point pattern by using clarkevans.test() of statspat.

ChildCare Services

The test hypotheses are:

Ho = The distribution of childcare services are randomly distributed.

H1= The distribution of childcare services are not randomly distributed. clustered

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("two.sided"),
                nsim=99)
## 
##  Clark-Evans test
##  No edge correction
##  Monte Carlo test based on 99 simulations of CSR with fixed n
## 
## data:  childcareSG_ppp
## R = 0.54015, p-value = 0.02
## alternative hypothesis: two-sided

Based on the P values(0.02) it is statistically significant and reject the Null hypotheses which means that the childcare services are not randomly distributed.

Based on the clark-evans test with 99 simulations, we are able to reject the null hypothesis at a 95% confidence level with a p-value of 0.02. , With the R-value of 0.53969 which is less than 1 shows that the distribution type of childcare services are clustered.

Kindergarten Services

The test hypotheses are:

Ho = The distribution of Kindergarten services are randomly distributed.

H1= The distribution of Kindergarten services are not randomly distributed.

clarkevans.test(KindergartenSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("two.sided"),
                nsim=99)
## 
##  Clark-Evans test
##  No edge correction
##  Monte Carlo test based on 99 simulations of CSR with fixed n
## 
## data:  KindergartenSG_ppp
## R = 0.50123, p-value = 0.02
## alternative hypothesis: two-sided

Based on the clark-evans test result (99 times simulation). P values(0.02) it is statistically significant and reject the Null hypotheses which means that the Kindergarten services are not randomly distributed.

Based on the clark evans test with 99 simulations, we are able to reject the null hypothesis at a 95% confidence level with a p-value of 0.02. , With the R-value of 0.50538 which is less than 1 shows that the distribution type of Kindergarten services are clustered.

2nd order spatial point patterns analysis

In this analysis, G function and L function will be used as the 2nd order spatial point patterns analysis.

G function indicate how does the services spaced in a point pattern.

Child Care Services

G_ChildCare = Gest(childcareSG_ppp, correction = "border")
plot(G_ChildCare)

From the G function plot above, it shows that increases rapidly at short distance, this is also an indication of clustering.Thus, this means that childcare services are spaced in a clustered point pattern.

Monte Carlo test with G-function to Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services in Singapore are randomly distributed.

H1= The distribution of childcare services in Singapore are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Confidence level is at 95%, and 1000 simulations are performed for this test, thus,t he alpha value is 0.001.

G_childcareRT<- envelope(childcareSG_ppp, Gest, nsim = 999)
## Generating 999 simulations of CSR  ...
## 1, 2, 3, ......10 [etd 48:49] .........20 [etd 48:53] .........
## 30 [etd 46:11] .........40 [etd 48:46] .........50 [etd 48:00] ........
## .60 [etd 46:32] .........70 [etd 45:10] .........80 [etd 44:12] .......
## ..90 [etd 44:32] .........100 [etd 43:44] .........110 [etd 42:50] ......
## ...120 [etd 41:58] .........130 [etd 41:10] .........140 [etd 40:20] .....
## ....150 [etd 39:40] .........160 [etd 38:55] .........170 [etd 38:16] ....
## .....180 [etd 37:44] .........190 [etd 37:09] .........200 [etd 36:37] ...
## ......210 [etd 36:04] .........220 [etd 35:31] .........230 [etd 35:00] ..
## .......240 [etd 34:28] .........250 [etd 33:56] .........260 [etd 33:28] .
## ........270 [etd 33:08] .........280 [etd 32:41] .........290
##  [etd 32:13] .........300 [etd 31:44] .........310 [etd 31:14] .........
## 320 [etd 30:48] .........330 [etd 30:19] .........340 [etd 29:49] ........
## .350 [etd 29:18] .........360 [etd 28:49] .........370 [etd 28:19] .......
## ..380 [etd 27:49] .........390 [etd 27:21] .........400 [etd 26:52] ......
## ...410 [etd 26:23] .........420 [etd 25:55] .........430 [etd 25:27] .....
## ....440 [etd 25:00] .........450 [etd 24:38] .........460 [etd 24:09] ....
## .....470 [etd 23:41] .........480 [etd 23:12] .........490 [etd 22:44] ...
## ......500 [etd 22:16] .........510 [etd 21:48] .........520 [etd 21:22] ..
## .......530 [etd 20:55] .........540 [etd 20:29] .........550 [etd 20:01] .
## ........560 [etd 19:33] .........570 [etd 19:06] .........580
##  [etd 18:45] .........590 [etd 18:20] .........600 [etd 17:52] .........
## 610 [etd 17:25] .........620 [etd 16:57] .........630 [etd 16:30] ........
## .640 [etd 16:02] .........650 [etd 15:35] .........660 [etd 15:08] .......
## ..670 [etd 14:41] .........680 [etd 14:14] .........690 [etd 13:47] ......
## ...700 [etd 13:24] .........710 [etd 13:03] .........720 [etd 12:39] .....
## ....730 [etd 12:14] .........740 [etd 11:49] .........750 [etd 11:24] ....
## .....760 [etd 10:58] .........770 [etd 10:32] .........780 [etd 10:06] ...
## ......790 [etd 9:39] .........800 [etd 9:12] .........810 [etd 8:46] ..
## .......820 [etd 8:19] .........830 [etd 7:52] .........840 [etd 7:25] .
## ........850 [etd 6:58] .........860 [etd 6:31] .........870
##  [etd 6:03] .........880 [etd 5:35] .........890 [etd 5:08] .........
## 900 [etd 4:40] .........910 [etd 4:12] .........920 [etd 3:44] ........
## .930 [etd 3:16] .........940 [etd 2:48] .........950 [etd 2:19] .......
## ..960 [etd 1:51] .........970 [etd 1:23] .........980 [etd 54 sec] ......
## ...990 [etd 26 sec] ........ 999.
## 
## Done.
plot(G_childcareRT)

With the use of Monte Carlo test with G-function on 999 times of simulation.We can concluded on the following

  • Estimated G(r) line is above the envelope this shows that it is statistically significant and so reject the null hypothesis.Thus, the childcare service are not randomly distributed.

  • The Estimated G(r) line is above the gray area (envelope) also indicate that the childcare services are in clustered pattern.

Kinder Services

G_kindergarten = Gest(KindergartenSG_ppp, correction = "border")
plot(G_kindergarten)

From the G function plot above, it shows that increases rapidly at short distance, this is also an indication of clustering.Thus, this means that kindergarten are spaced in a clustered point pattern.

Monte Carlo test with G-function to Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of Kindergarten services in Singapore are randomly distributed.

H1= The distribution of Kindergarten services in Singapore are not randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001.

Confidence level is at 95%, and 1000 simulations are performed for this test, thus,t he alpha value is 0.001.

G_kindergartenRT<- envelope(KindergartenSG_ppp, Gest, nsim = 999)
## Generating 999 simulations of CSR  ...
## 1, 2, 3, ......10 [etd 23:41] .........20 [etd 21:48] .........
## 30 [etd 21:33] .........40 [etd 20:53] .........50 [etd 20:28] ........
## .60 [etd 21:30] .........70 [etd 22:03] .........80 [etd 22:27] .......
## ..90 [etd 22:38] .........100 [etd 22:45] .........110 [etd 22:50] ......
## ...120 [etd 22:26] .........130 [etd 22:08] .........140 [etd 22:05] .....
## ....150 [etd 21:59] .........160 [etd 21:55] .........170 [etd 21:47] ....
## .....180 [etd 21:38] .........190 [etd 21:30] .........200 [etd 21:21] ...
## ......210 [etd 20:44] .........220 [etd 20:09] .........230 [etd 19:35] ..
## .......240 [etd 19:05] .........250 [etd 18:37] .........260 [etd 18:10] .
## ........270 [etd 17:44] .........280 [etd 17:21] .........290
##  [etd 16:57] .........300 [etd 16:34] .........310 [etd 16:11] .........
## 320 [etd 15:51] .........330 [etd 15:34] .........340 [etd 15:15] ........
## .350 [etd 14:59] .........360 [etd 14:42] .........370 [etd 14:25] .......
## ..380 [etd 14:09] .........390 [etd 13:52] .........400 [etd 13:36] ......
## ...410 [etd 13:20] .........420 [etd 13:04] .........430 [etd 12:49] .....
## ....440 [etd 12:34] .........450 [etd 12:19] .........460 [etd 12:04] ....
## .....470 [etd 11:49] .........480 [etd 11:35] .........490 [etd 11:20] ...
## ......500 [etd 11:06] .........510 [etd 10:52] .........520 [etd 10:38] ..
## .......530 [etd 10:24] .........540 [etd 10:10] .........550 [etd 9:56] .
## ........560 [etd 9:42] .........570 [etd 9:28] .........580
##  [etd 9:14] .........590 [etd 9:00] .........600 [etd 8:46] .........
## 610 [etd 8:32] .........620 [etd 8:18] .........630 [etd 8:04] ........
## .640 [etd 7:51] .........650 [etd 7:37] .........660 [etd 7:24] .......
## ..670 [etd 7:10] .........680 [etd 6:57] .........690 [etd 6:43] ......
## ...700 [etd 6:30] .........710 [etd 6:17] .........720 [etd 6:03] .....
## ....730 [etd 5:49] .........740 [etd 5:35] .........750 [etd 5:22] ....
## .....760 [etd 5:08] .........770 [etd 4:54] .........780 [etd 4:41] ...
## ......790 [etd 4:28] .........800 [etd 4:14] .........810 [etd 4:01] ..
## .......820 [etd 3:48] .........830 [etd 3:34] .........840 [etd 3:21] .
## ........850 [etd 3:08] .........860 [etd 2:55] .........870
##  [etd 2:42] .........880 [etd 2:30] .........890 [etd 2:17] .........
## 900 [etd 2:04] .........910 [etd 1:51] .........920 [etd 1:38] ........
## .930 [etd 1:26] .........940 [etd 1:13] .........950 [etd 1:01] .......
## ..960 [etd 48 sec] .........970 [etd 36 sec] .........980 [etd 23 sec] ......
## ...990 [etd 11 sec] ........ 999.
## 
## Done.
plot(G_kindergartenRT)

With the use of Monte Carlo test with G-function on 999 times of simulation.We can concluded on the following

  • Estimated G(r) line is above the envelope this shows that it is statistically significant and so reject the null hypothesis.Thus, the kindergarten service are not randomly distributed.

  • The Estimated G(r) line is above the gray area (envelope) also indicate that the kindergarten services are in clustered pattern.

Thus based on the above studies, we can conclude that the point pattern for both child care services and kindergarten services are clustered and not randomly distributed around Singapore. Which it fit with our daily live observation, as child care and kindergarten are frequently build around the residential area.

Kernel Density Map

In this section, I will be Computing kernel density by using density() function.

Rescale

childcareSG_ppp.km <- rescale(childcareSG_ppp, 1000, "km")
KindergartenSG_ppp.km <- rescale(KindergartenSG_ppp, 1000, "km")

ChildCare KDE

Based on the above statistical analysis both childcare services and kindergarten service are tight clusters distribution.Hence bw.ppl() algorithm is used in this analysis as it produce the more appropriate values when the pattern consists predominantly of tight clusters. The reason bw.diggle() function is not used in this analysis is because the diggle function is to detect a single tight cluster in the midst of random noise. So based on the distribution point that show in the Combining service point with SG Polygon section shows that the cluster is not a single tight cluster in the midst of random noise. Hence, bw.ppl() function is selected.

kde_childcareSG.ppl <- density(childcareSG_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="gaussian")

bw.ppl(childcareSG_ppp.km)
##     sigma 
## 0.3072432
plot(kde_childcareSG.ppl, main = "bw.ppl")

Based on the Childcare KDE map, child care services are clustered within the same region and not randomly distributed.

Kindergarten Service KDE

kde_KindergartenSG.ppl <- density(KindergartenSG_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="gaussian")

bw.ppl(KindergartenSG_ppp.km)
##    sigma 
## 0.611329
plot(kde_KindergartenSG.ppl, main = "bw.ppl")

Based on the Childcare KDE map, kindergarten services are clustered within the same region and not randomly distributed.

OSM with Raster Layer

In order to plot the KDE map on to the Open Street Map, it is required to perform the following step

    1. convert KDE output to Spatial Grid Data frame ojbect.
    1. Convert Spatial Grid Data frame object into raster object.
    1. Assigning the projection system.

Converting KDE output –> Spatial Grid Dataframe

gridded_kde_childcareSG <- as.SpatialGridDataFrame.im(kde_childcareSG.ppl)

gridded_kde_KindergartenSG <- as.SpatialGridDataFrame.im(kde_KindergartenSG.ppl)

Converting gridded kernal density objects –> RasterLayer object.

gridded_kde_childcareSG_raster <- raster(gridded_kde_childcareSG )
gridded_kde_KindergartenSG_raster <- raster(gridded_kde_KindergartenSG)

gridded_kde_childcareSG_raster
## class      : RasterLayer 
## dimensions : 128, 128, 16384  (nrow, ncol, ncell)
## resolution : 0.419757, 0.2695907  (x, y)
## extent     : 2.667538, 56.39644, 15.74872, 50.25633  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : v 
## values     : -1.055861e-14, 32.52153  (min, max)
gridded_kde_KindergartenSG_raster
## class      : RasterLayer 
## dimensions : 128, 128, 16384  (nrow, ncol, ncell)
## resolution : 0.419757, 0.2695907  (x, y)
## extent     : 2.667538, 56.39644, 15.74872, 50.25633  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : v 
## values     : -1.507226e-15, 6.372976  (min, max)

Assigning project systems

projection(gridded_kde_childcareSG_raster) <- CRS("+init=EPSG:3414 +datum=WGS84 +units=km")
projection(gridded_kde_KindergartenSG_raster) <- CRS("+init=EPSG:3414 +datum=WGS84 +units=km")
gridded_kde_childcareSG_raster
## class      : RasterLayer 
## dimensions : 128, 128, 16384  (nrow, ncol, ncell)
## resolution : 0.419757, 0.2695907  (x, y)
## extent     : 2.667538, 56.39644, 15.74872, 50.25633  (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 +datum=WGS84 +units=km +no_defs 
## source     : memory
## names      : v 
## values     : -1.055861e-14, 32.52153  (min, max)
gridded_kde_KindergartenSG_raster
## class      : RasterLayer 
## dimensions : 128, 128, 16384  (nrow, ncol, ncell)
## resolution : 0.419757, 0.2695907  (x, y)
## extent     : 2.667538, 56.39644, 15.74872, 50.25633  (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 +datum=WGS84 +units=km +no_defs 
## source     : memory
## names      : v 
## values     : -1.507226e-15, 6.372976  (min, max)

Open Stree Map with Raster Layer

This section will be combining the open street map with the following 2 raster layer:

  • Child Care Services
  • Kindergarten Services

Child Care Services & Kindergarten Services

tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap")+
tm_shape(gridded_kde_childcareSG_raster) + 
  tm_raster("v", n=5, title= "Child Care Services", alpha=0.4 ,drop.levels = TRUE, palette = "Reds")+ 
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE) +
tm_shape(gridded_kde_KindergartenSG_raster) + 
  tm_raster("v", n=5, title= "Kindergarten Services",breaks=c(0,5,10,15,20,25,30,35), alpha=0.2,drop.levels = TRUE, palette = "Reds")
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Warning: Values have found that are less than the lowest break

Based on the KDE on osm map, we can conclude that both childcare and kindergarten services are highly correlated as both services located at the same region.However, the number of services provided and density level for child care services are higher than kindergarten services.

On the OSM with the raster layer, also shows that the higher density for both services are located at residential area such as: SengKang, Bungkok,Serangoon,woodland and Sambawang. This also reflected based on what is happening in current. As those area mentioned above is either already matured town or going to be a matured town where there are a lot of BTO building is going on or its a young town where young family is located.

Compare the advantages of kernel density maps over point maps.

  • The Kernel density map shows clearer on the intensity of the area, where as it is not so easily to observed over the point maps.

  • With the Density based techniques it characterize the pattern in terms of its distribution vis-a-vis the study area–a first-order property of the pattern.