pacman::p_load(sf,tidyverse)Hands-on Exercise 01 Data Wrangling with R
Overview
Hand-on exercise 1 teaches how to import and wrangle geospatial data using R packages.
Getting Started
Importing Geospatial Data
Import polygon feature data
mpsz <- st_read(dsn="geospatial",layer="MP14_SUBZONE_WEB_PL")Reading layer `MP14_SUBZONE_WEB_PL' from data source
`C:\Soe Htet\ISSS624\Hands-on Exercise 01\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21cyclingpath <- st_read(dsn="geospatial",
layer="CyclingPath")Reading layer `CyclingPath' from data source
`C:\Soe Htet\ISSS624\Hands-on Exercise 01\geospatial' using driver `ESRI Shapefile'
Simple feature collection with 1625 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 12711.19 ymin: 28711.33 xmax: 42626.09 ymax: 48948.15
Projected CRS: SVY21preschool <- st_read("geospatial/pre-schools-location-kml.kml")Reading layer `PRESCHOOLS_LOCATION' from data source
`C:\Soe Htet\ISSS624\Hands-on Exercise 01\geospatial\pre-schools-location-kml.kml'
using driver `KML'
Simple feature collection with 1359 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84Working with st_geometry()
st_geometry(mpsz)Geometry set for 323 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
First 5 geometries:MULTIPOLYGON (((31495.56 30140.01, 31980.96 296...MULTIPOLYGON (((29092.28 30021.89, 29119.64 300...MULTIPOLYGON (((29932.33 29879.12, 29947.32 298...MULTIPOLYGON (((27131.28 30059.73, 27088.33 297...MULTIPOLYGON (((26451.03 30396.46, 26440.47 303...Working with glimpse()
glimpse(mpsz)Rows: 323
Columns: 16
$ OBJECTID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ SUBZONE_NO <int> 1, 1, 3, 8, 3, 7, 9, 2, 13, 7, 12, 6, 1, 5, 1, 1, 3, 2, 2, …
$ SUBZONE_N <chr> "MARINA SOUTH", "PEARL'S HILL", "BOAT QUAY", "HENDERSON HIL…
$ SUBZONE_C <chr> "MSSZ01", "OTSZ01", "SRSZ03", "BMSZ08", "BMSZ03", "BMSZ07",…
$ CA_IND <chr> "Y", "Y", "Y", "N", "N", "N", "N", "Y", "N", "N", "N", "N",…
$ PLN_AREA_N <chr> "MARINA SOUTH", "OUTRAM", "SINGAPORE RIVER", "BUKIT MERAH",…
$ PLN_AREA_C <chr> "MS", "OT", "SR", "BM", "BM", "BM", "BM", "SR", "QT", "QT",…
$ REGION_N <chr> "CENTRAL REGION", "CENTRAL REGION", "CENTRAL REGION", "CENT…
$ REGION_C <chr> "CR", "CR", "CR", "CR", "CR", "CR", "CR", "CR", "CR", "CR",…
$ INC_CRC <chr> "5ED7EB253F99252E", "8C7149B9EB32EEFC", "C35FEFF02B13E0E5",…
$ FMEL_UPD_D <date> 2014-12-05, 2014-12-05, 2014-12-05, 2014-12-05, 2014-12-05…
$ X_ADDR <dbl> 31595.84, 28679.06, 29654.96, 26782.83, 26201.96, 25358.82,…
$ Y_ADDR <dbl> 29220.19, 29782.05, 29974.66, 29933.77, 30005.70, 29991.38,…
$ SHAPE_Leng <dbl> 5267.381, 3506.107, 1740.926, 3313.625, 2825.594, 4428.913,…
$ SHAPE_Area <dbl> 1630379.27, 559816.25, 160807.50, 595428.89, 387429.44, 103…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((31495.56 30..., MULTIPOLYGON (…Working with head()
head(mpsz,n=5)Simple feature collection with 5 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 25867.68 ymin: 28369.47 xmax: 32362.39 ymax: 30435.54
Projected CRS: SVY21
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
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
Y_ADDR SHAPE_Leng SHAPE_Area geometry
1 29220.19 5267.381 1630379.3 MULTIPOLYGON (((31495.56 30...
2 29782.05 3506.107 559816.2 MULTIPOLYGON (((29092.28 30...
3 29974.66 1740.926 160807.5 MULTIPOLYGON (((29932.33 29...
4 29933.77 3313.625 595428.9 MULTIPOLYGON (((27131.28 30...
5 30005.70 2825.594 387429.4 MULTIPOLYGON (((26451.03 30...Plotting the geospatial data
plot(mpsz)Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all
plot(st_geometry(mpsz))
plot(mpsz["PLN_AREA_N"])
Working with Projection
One of the common issue that can happen during importing geospatial data into R is that the coordinate system of the source data was either missing (such as due to missing .proj for ESRI shapefile) or wrongly assigned during the importing process.
This is an example the coordinate system of mpsz simple feature data frame by using st_crs() of sf package as shown in the code chunk below.
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]]]]mpsz3414 = st_set_crs(mpsz,3414)Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
thatst_crs(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["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]Transforming the projection of preschool from wgs84 to svy21
preschool3414 = st_transform(preschool,crs=3414)
preschool3414Simple feature collection with 1359 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 11203.01 ymin: 25667.6 xmax: 45404.24 ymax: 49300.88
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name
1 kml_1
2 kml_2
3 kml_3
4 kml_4
5 kml_5
6 kml_6
7 kml_7
8 kml_8
9 kml_9
10 kml_10
Description
1 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>BIG FOOT PRE SCHOOL LLP</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT9281</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>196, WEST COAST ROAD, SINGAPORE 127375</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>127375</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>838CD358794FD031</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
2 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>POSSO PRESCHOOL @ WEST COAST RISE PTE LTD</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT8684</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>30, WEST COAST RISE, HONG LEONG GARDEN, SINGAPORE 127473</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>127473</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>F331CEB175F9C254</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
3 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>GENESIS CHILD CARE PTE. LTD.</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT9132</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>2A, JUBILEE ROAD, SINGAPORE 128524</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>128524</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4C2E7E55019A633F</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
4 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>LITTLE FOOTPRINTS PRESCHOOL PTE. LTD.</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT9260</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>6, JUBILEE ROAD, SINGAPORE 128531</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>128531</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>DDF98422A198387B</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
5 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>AMAR KIDZ @ WEST COAST LLP</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT9016</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>8, JALAN LEMPENG, #02 - 03, PARK WEST CONDO, SINGAPORE 128796</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>128796</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>EAB3263D23F126AF</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
6 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>TCC PRESCHOOL FABER PTE LTD</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT9299</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>60, FABER TERRACE, FABER HILLS, SINGAPORE 129040</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>129040</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>195E3739B77E6A5F</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
7 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>ACEKIDZ @ COMMUNITY</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>PT5950</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>2, CLEMENTI WEST ST 2, #03 - 06, WEST COAST COMMUNITY CENTRE, SINGAPORE 129605</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>129605</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>9B1070EE1CB4A3E2</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
8 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>PCF SPARKLETOTS PRESCHOOL @ QUEENSTOWN BLK 145 (CC)</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>ST0092</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>145, MEI LING STREET, #01 - 137, SINGAPORE 140145</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>140145</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>820E90716985CCCA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
9 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>PCF SPARKLETOTS PRESCHOOL @ QUEENSTOWN BLK 53A (CC)</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>ST0176</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>53A, STRATHMORE AVENUE, #01 - 01, FORFAR HEIGHTS, SINGAPORE 143053</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>143053</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>A7DC7D2C961A8822</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
10 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>CENTRE_NAME</th> <td>MY FIRST SKOOL</td> </tr><tr bgcolor=""> <th>CENTRE_CODE</th> <td>NT0510</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESS</th> <td>106, HENDERSON CRESCENT, #01 - 37, SINGAPORE 150106</td> </tr><tr bgcolor=""> <th>POSTAL_CODE</th> <td>150106</td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>EB3942B460BB5CBC</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20171208174048</td> </tr></table></center>
geometry
1 POINT Z (19997.26 32333.17 0)
2 POINT Z (19126.75 33114.35 0)
3 POINT Z (20345.12 31934.56 0)
4 POINT Z (20400.31 31952.36 0)
5 POINT Z (19810.78 33140.31 0)
6 POINT Z (19550.92 33770.18 0)
7 POINT Z (20378.07 31665.55 0)
8 POINT Z (24835.77 30689.38 0)
9 POINT Z (25139.3 30636.01 0)
10 POINT Z (26771.14 30203.71 0)Importing and Converting An Aspatial Data
listings = read.csv("aspatial/listings.csv")listings_sf <- st_as_sf (listings,
coords=c("longitude","latitude"),
crs = 4326) %>%
st_transform(crs=3414)glimpse(listings_sf)Rows: 4,252
Columns: 15
$ id <int> 50646, 71609, 71896, 71903, 275343, 275…
$ name <chr> "Pleasant Room along Bukit Timah", "Ens…
$ host_id <int> 227796, 367042, 367042, 367042, 1439258…
$ host_name <chr> "Sujatha", "Belinda", "Belinda", "Belin…
$ neighbourhood_group <chr> "Central Region", "East Region", "East …
$ neighbourhood <chr> "Bukit Timah", "Tampines", "Tampines", …
$ room_type <chr> "Private room", "Private room", "Privat…
$ price <int> 80, 178, 81, 81, 52, 40, 72, 41, 49, 49…
$ minimum_nights <int> 90, 90, 90, 90, 14, 14, 90, 8, 14, 14, …
$ number_of_reviews <int> 18, 20, 24, 48, 20, 13, 133, 105, 14, 1…
$ last_review <chr> "2014-07-08", "2019-12-28", "2014-12-10…
$ reviews_per_month <dbl> 0.22, 0.28, 0.33, 0.67, 0.20, 0.16, 1.2…
$ calculated_host_listings_count <int> 1, 4, 4, 4, 50, 50, 7, 1, 50, 50, 50, 4…
$ availability_365 <int> 365, 365, 365, 365, 353, 364, 365, 90, …
$ geometry <POINT [m]> POINT (22646.02 35167.9), POINT (…Geoprocessing with sf package
Besides providing functions to handling (i.e. importing, exporting, assigning projection, transforming projection etc) geospatial data, sf package also offers a wide range of geoprocessing (also known as GIS analysis) functions.
In this section, you will learn how to perform two commonly used geoprocessing functions, namely buffering and point in polygon count.
Buffering
The scenario:
The authority is planning to upgrade the exiting cycling path. To do so, they need to acquire 5 metres of reserved land on the both sides of the current cycling path. You are tasked to determine the extend of the land need to be acquired and their total area.
The solution:
Firstly, st_buffer() of sf package is used to compute the 5-meter buffers around cycling paths
buffer_cycling <- st_buffer(cyclingpath,dist=5,nQuadSegs = 30)buffer_cycling$AREA <- st_area(buffer_cycling)sum(buffer_cycling$AREA)773143.9 [m^2]Point in Polygon Count
The scenario:
A pre-school service group want to find out the numbers of pre-schools in each Planning Subzone.
The solution:
The code chunk below performs two operations at one go. Firstly, identify pre-schools located inside each Planning Subzone by using st_intersects(). Next, length() of Base R is used to calculate numbers of pre-schools that fall inside each planning subzone.
mpsz3414$`PreSch Count`<- lengths(st_intersects(mpsz3414, preschool3414))summary(mpsz3414$`PreSch Count`) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 2.000 4.207 6.000 37.000 top_n(mpsz3414, 1, `PreSch Count`)Simple feature collection with 1 feature and 16 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 23449.05 ymin: 46001.23 xmax: 25594.22 ymax: 47996.47
Projected CRS: SVY21 / Singapore TM
OBJECTID SUBZONE_NO SUBZONE_N SUBZONE_C CA_IND PLN_AREA_N PLN_AREA_C
1 290 3 WOODLANDS EAST WDSZ03 N WOODLANDS WD
REGION_N REGION_C INC_CRC FMEL_UPD_D X_ADDR Y_ADDR
1 NORTH REGION NR C90769E43EE6B0F2 2014-12-05 24506.64 46991.63
SHAPE_Leng SHAPE_Area geometry PreSch Count
1 6603.608 2553464 MULTIPOLYGON (((24786.75 46... 37mpsz3414$Area <- mpsz3414 %>%
st_area()mpsz3414 <- mpsz3414 %>%
mutate(`PreSch Density` = `PreSch Count`/Area * 1000000)Exploratory Data Analysis
hist(mpsz3414$`PreSch Density`)
With ggplot2
ggplot(data=mpsz3414,aes(x=as.numeric(`PreSch Density`)))+
geom_histogram(bins=20,color='black',fill='light blue')+
labs(title = "Are pre-school evenly distributed in Singapore?",
subtitle='...',
x="Pre sch density (per km sq)",
y="Frequency")
ggplot(data=mpsz3414,aes(x=as.numeric(`PreSch Density`),y=`PreSch Count`))+
geom_point()+
labs(title = "Are pre-school evenly distributed in Singapore?",
subtitle='...',
x="Pre sch density (per km sq)",
y="Count")
2 Choropleth Mapping with R
pacman::p_load(tmap)Importing geospatial data into R
popdata <- read.csv('aspatial/respopagesextod2011to2020.csv')Data Wrangling
popdata2020 <- popdata %>%
filter(Time == 2020) %>%
group_by(PA, SZ, AG) %>%
summarise(POP = sum(Pop)) %>%
ungroup()%>%
pivot_wider(names_from=AG,
values_from=POP)%>%
mutate(YOUNG = rowSums(.[3:6])
+rowSums(.[12])) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[7:11])+
rowSums(.[13:15]))%>%
mutate(AGED=rowSums(.[16:21])) %>%
mutate(TOTAL=rowSums(.[3:21])) %>%
mutate(DEPENDENCY = (`YOUNG` + `AGED`)/`ECONOMY ACTIVE`) %>%
select(PA, SZ, YOUNG,
`ECONOMY ACTIVE`, AGED,
TOTAL, DEPENDENCY)`summarise()` has grouped output by 'PA', 'SZ'. You can override using the
`.groups` argument.Joining the attribute data and geospatial data
popdata2020 <- popdata2020 %>%
mutate_at(.vars=vars(PA,SZ),
.funs=funs(toupper)) %>%
filter (`ECONOMY ACTIVE`>0)Warning: `funs()` was deprecated in 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_lifecycle_warnings()` to see where this warning was generated.mpsz_pop2020 <- left_join(mpsz,popdata2020,by=c("SUBZONE_N" = "SZ"))write_rds(mpsz_pop2020, "rds/mpszpop2020.rds")Choropleth Mapping Geospatial Data Using qtm()
Two approaches can be used to prepare thematic map using tmap, they are:
Plotting a thematic map quickly by using qtm().
Plotting highly customisable thematic map by using tmap elements.
Plotting a choropleth map using qtm()
library(tmap)
tmap_mode("plot")tmap mode set to plottingqtm(mpsz_pop2020,
fill = "DEPENDENCY")
Creating a choropleth map by using tmap's elements
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
title = "Dependency ratio") +
tm_layout(main.title = "Distribution of Dependency Ratio by subzone",
main.title.position = "center",
main.title.size = 1.2,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 2) +
tm_scale_bar() +
tm_grid(alpha =0.2)
Drawing a base map
tm_shape(mpsz_pop2020) +
tm_polygons()
Drawing a choropleth map using tm_polygons()
tm_shape(mpsz_pop2020)+
tm_polygons("DEPENDENCY")
Drawing a choropleth map using tm_fill() & tm_borders()
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY")
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY")+
tm_borders(lwd = 0.1, alpha = 1)
Data Classification Methods of tmap
Plotting choropleth maps with built in classification methods
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
n = 5,
style = "jenks") +
tm_borders(alpha = 0.5)
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
n = 5,
style = "equal") +
tm_borders(alpha = 0.5)
Plotting Choropleth Maps With Custom Breaks
summary(mpsz_pop2020$DEPENDENCY) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0000 0.6519 0.7025 0.7742 0.7645 19.0000 92 tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
breaks = c(0, 0.60, 0.70, 0.80, 0.90, 1.00)) +
tm_borders(alpha = 0.5)Warning: Values have found that are higher than the highest break
Colour Schemes
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "-Greens") +
tm_borders(alpha = 0.5)
Map Layout
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
style = "jenks",
palette = "Blues",
legend.hist = TRUE,
legend.is.portrait = TRUE,
legend.hist.z = 0.1) +
tm_layout(main.title = "Distribution of Dependency Ratio by planning subzone \n(Jenks classification)",
main.title.position = "center",
main.title.size = 1,
legend.height = 0.45,
legend.width = 0.35,
legend.outside = FALSE,
legend.position = c("right", "bottom"),
frame = FALSE) +
tm_borders(alpha = 0.5)
Map Style
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "-Greens") +
tm_borders(alpha = 0.5) +
tmap_style("classic")tmap style set to "classic"other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor" 
Cartographic furniture
tm_shape(mpsz_pop2020)+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
title = "No. of persons") +
tm_layout(main.title = "Distribution of Dependency Ratio \nby planning subzone",
main.title.position = "center",
main.title.size = 1.2,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 2) +
tm_scale_bar(width = 0.15) +
tm_grid(lwd = 0.1, alpha = 0.2) +
tm_credits("Source: Planning Sub-zone boundary from Urban Redevelopment Authorithy (URA)\n and Population data from Department of Statistics DOS",
position = c("left", "bottom"))
Drawing Small Multiple Choropleth Maps
tm_shape(mpsz_pop2020)+
tm_fill(c("YOUNG", "AGED"),
style = "equal",
palette = "Blues") +
tm_layout(legend.position = c("right", "bottom")) +
tm_borders(alpha = 0.5) +
tmap_style("white")tmap style set to "white"other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor" 
tm_shape(mpsz_pop2020)+
tm_polygons(c("DEPENDENCY","AGED"),
style = c("equal", "quantile"),
palette = list("Blues","Greens")) +
tm_layout(legend.position = c("right", "bottom"))
youngmap <- tm_shape(mpsz_pop2020)+
tm_polygons("YOUNG",
style = "quantile",
palette = "Blues")
agedmap <- tm_shape(mpsz_pop2020)+
tm_polygons("AGED",
style = "quantile",
palette = "Blues")
tmap_arrange(youngmap, agedmap, asp=1, ncol=2)
Mapping spatial object meeting a certain criterion
tm_shape(mpsz_pop2020[mpsz_pop2020$REGION_N=="CENTRAL REGION", ])+
tm_fill("DEPENDENCY",
style = "quantile",
palette = "Blues",
legend.hist = TRUE,
legend.is.portrait = TRUE,
legend.hist.z = 0.1) +
tm_layout(legend.outside = TRUE,
legend.height = 0.45,
legend.width = 5.0,
legend.position = c("right", "bottom"),
frame = FALSE) +
tm_borders(alpha = 0.5)Warning in pre_process_gt(x, interactive = interactive, orig_crs =
gm$shape.orig_crs): legend.width controls the width of the legend within a map.
Please use legend.outside.size to control the width of the outside legend