Geospatial analysis
Importing the data
library(terra)
library(tidyterra)
library(sf)
library(tidyverse)
library(patchwork)
mrtsg.cs <- read.csv("mrtsg.csv") # importing our mrt data
mrt <- st_as_sf(mrtsg.cs, coords = c("Longitude", "Latitude"), crs = 4326) # converting the data to sf
# importing the shapefiles and raster
res_shp <- st_read("data (2)/data/master_plan/MP2019zoning_residential_buffer.shp")
## Reading layer `MP2019zoning_residential_buffer' from data source
## `C:\Users\mikek\OneDrive - University College Dublin\Documents\Tutorial 4\data (2)\data\master_plan\MP2019zoning_residential_buffer.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 103.6826 ymin: 1.237671 xmax: 103.9898 ymax: 1.462815
## Geodetic CRS: WGS 84
sgp_boundary <- st_read("data (2)/data/gadm36_SGP_shp/gadm36_SGP_0.shp")
## Reading layer `gadm36_SGP_0' from data source
## `C:\Users\mikek\OneDrive - University College Dublin\Documents\Tutorial 4\data (2)\data\gadm36_SGP_shp\gadm36_SGP_0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 103.6091 ymin: 1.16639 xmax: 104.0858 ymax: 1.471388
## Geodetic CRS: WGS 84
rail.data <- st_read("data (2)/data/osm/osm_singapore_rail.shp")
## Reading layer `osm_singapore_rail' from data source
## `C:\Users\mikek\OneDrive - University College Dublin\Documents\Tutorial 4\data (2)\data\osm\osm_singapore_rail.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1320 features and 7 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 103.6291 ymin: 1.250935 xmax: 103.9919 ymax: 1.452361
## Geodetic CRS: WGS 84
regions <- st_read("data (2)/data/planning-area-census2010-shp/Planning_Area_Census2010.shp")
## Reading layer `Planning_Area_Census2010' from data source
## `C:\Users\mikek\OneDrive - University College Dublin\Documents\Tutorial 4\data (2)\data\planning-area-census2010-shp\Planning_Area_Census2010.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 55 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2637.869 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## Projected CRS: SVY21
map_500m <- rast("data (2)/data/SGP_OLI8_500m.tif")
map of all layers together
# plotting all shapefiles on the raster
ggplot() +
geom_spatraster_rgb(data = map_500m - 5000, r = 4, g = 5, b = 6, max_col_value = 12000) +
geom_sf(data = res_shp, aes(color = 'Residential areas'), fill = NA, size = 0.25) +
geom_sf(data = rail.data, aes(color = "MRT lines")) +
geom_sf(data = mrt, aes(color = "MRT stations"), shape = "x", size = 1.5) +
scale_color_manual("",values = c("purple","aquamarine","white")) +
theme(legend.position = c(0.8,0.3), legend.background = element_blank(),
legend.text = element_text(color = "white"), legend.key = element_blank())
Difference in greenness and temperature ny area
# making our data usable for the map
map_2km <- aggregate(map_500m, 4) # aggregating to 2km pixels
map_2km$NDVI <- (map_2km$NIR - map_2km$Red) / (map_2km$NIR + map_2km$Red) # adding NDVI band
regions <- regions %>% st_transform(crs = 4326)
data.census <- rasterize(vect(regions), map_2km, field = "PLN_AREA_N") # using rasterize to convert layers to pixels
names(map_2km) # making sure column names are valid
## [1] "Sfc_Temp_C" "SWIR" "NIR" "Red" "Green"
## [6] "Blue" "NDVI"
cells <- as.data.frame(c(map_2km[[c("Sfc_Temp_C","NDVI")]], data.census)) # converting to a dataframe to get mean of each area
cell_area <- cells %>% group_by(PLN_AREA_N) %>%
summarize(Sfc_Temp_C = mean(Sfc_Temp_C, na.rm = TRUE),
NDVI = mean(NDVI, na.rm = TRUE)) # get mean area of each region
census <- regions %>% left_join(cell_area, by = "PLN_AREA_N") # rejoining with the shapefile
# mapping the data
base <- ggplot() +
coord_sf(xlim = c(103, 104), ylim = c(1.1, 1.4)) # setting the x and y limits of our map
b_green <- base +
scale_fill_viridis_c("NDVI", limits = c(-0.16, 0.5)) +
ggtitle("Greenness")
b_temp <- base +
scale_fill_viridis_c("Surface\nTemp (C)", option = "B", limits = c(25, 56)) +
ggtitle("Surface temperature")
temp.data <- b_temp +
geom_sf(data = census,aes(fill = Sfc_Temp_C), color = "white", size = 0.25)
green.data <- b_green +
geom_sf(data = census, aes(fill = NDVI), color = "white", size = 0.25)
m_temp_rast <- b_temp + geom_spatraster(data = map_2km, aes(fill = Sfc_Temp_C))
m_greenness_rast <- b_green + geom_spatraster(data = map_2km, aes(fill = NDVI)) # plots surface temperature and greenness to a map of singapore
b_green <- base +
scale_fill_viridis_c("NDVI", limits = c(-0.16, 0.5)) +
ggtitle("Greenness")
b_temp <- base +
scale_fill_viridis_c("Surface\nTemp (C)", option = "B", limits = c(25, 56)) +
ggtitle("Surface temperature")
temp.data <- b_temp +
geom_sf(data = census,aes(fill = Sfc_Temp_C), color = "white", size = 0.25)
green.data <- b_green +
geom_sf(data = census, aes(fill = NDVI), color = "white", size = 0.25)
m_temp_rast <- b_temp + geom_spatraster(data = map_2km, aes(fill = Sfc_Temp_C))
m_greenness_rast <- b_green + geom_spatraster(data = map_2km, aes(fill = NDVI))
relationship between temperature and greenness
map_px <- as.data.frame(map_2km, xy = TRUE) # using xy = TRUE so that are data has and x and y values
ggplot(map_px) + geom_point(aes(NDVI, Sfc_Temp_C)) # clearly a positive correlation in general between temp and greenness but some negative correlation between the middle to right
map_px <- map_px %>%
mutate(below_line = Sfc_Temp_C < NDVI * -55 + 52) # mutating to have separate data for values temp < green
p_b_line <- ggplot(map_px) +
geom_point(aes(NDVI, Sfc_Temp_C, color = below_line)) + ggtitle("Separate Line") # data to graph above and below the line
m_b_line <- ggplot(map_px) +
geom_raster(aes(x, y, fill = below_line)) + coord_equal() + ggtitle("Map of separate points") # data to map the separate lines
p_b_line
m_b_line # map shows that Singapore land is the main area where data is above the line
identifying water pixels
map_px <- map_px %>%
mutate(MNDWI = (Green - SWIR) / (Green + SWIR),
water = MNDWI > 0.0)
p_mndwi <- ggplot(map_px) + geom_point(aes(NDVI, Sfc_Temp_C, color = water), alpha = 0.2)
m_mndwi <- ggplot(map_px) +
geom_raster(aes(x, y, fill = MNDWI)) +
coord_equal()
p_mndwi # This shows that the areas that are water have lower greenery and temperatures. There are still some land points with the same values which is interesting
m_mndwi # map further strengthens this argument
showing urban area and park area temperature and green
variability
graph_green <- ggplot(map_px) + geom_point(aes(NDVI, Sfc_Temp_C, color = NDVI), alpha = 0.2) +
geom_point(data = map_px %>% filter(!water & below_line) %>% mutate(bands = "land_below_line"),
aes(NDVI, Sfc_Temp_C), color = "green")
map_green_b <- ggplot(map_px %>% filter(!water & below_line)) +
geom_raster(data = map_px, aes(x, y, fill = NDVI)) +
scale_fill_viridis_c()+
geom_raster(aes(x, y), fill = "red", alpha = 0.5) + coord_equal() +
ggtitle("NVDI & points below line") +
labs(caption = "Red indicates below the main cluster and classified as land (MNDWI < 0)")
graph_green +
map_green_b # can clearly see that more urban areas have higher temperatures and lower greenery meaning there is implications for reflections from buildings
Livability based on temperatures and greenery
# mapping the four measures
res_df <- rasterize(vect(res_shp), map_2km) # rasterizing our res data
r_data <- c(map_2km, res_df)
names(r_data) <- c(names(map_2km), "residential")
r_df_m <- as.data.frame(r_data, xy = TRUE) %>% as_tibble()
r_cell <- r_df_m %>% filter(!is.na(residential)) %>%
st_as_sf(coords = c("x","y"), crs = 4326) # filtering to just the res cells
distance <- st_distance(r_cell, mrt) %>% as_tibble()
r_cell$mrt_dist <- apply(distance, 1, min) # distance to the MRT while extracting minimums
scale_0_to_1 <- function(x, na.rm = TRUE) {
x_scaled <- (x - min(x, na.rm = na.rm)) / (max(x, na.rm = na.rm) - min(x, na.rm = na.rm))
return(x_scaled)
}
r_df_m2 <- r_df_m %>%
left_join(r_cell %>% st_set_geometry(NULL)) %>%
mutate(TemperatureIndex = 1 - scale_0_to_1(Sfc_Temp_C),
GreennessIndex = scale_0_to_1(NDVI),
TransportIndex = 1 - scale_0_to_1(mrt_dist)) %>%
dplyr::select(x, y, ends_with("Index")) # calculating the indices
df_long <- r_df_m2 %>% pivot_longer(ends_with("Index"), names_to = "Variables",values_to = "value") # pivoting òur data longer
df_live <- df_long %>%
group_by(x,y) %>%
summarize(variable = "Livability", value = mean(value)) # calculating livability index
index <- bind_rows(df_long, df_live) # join living index with indices
ggplot(index) +
geom_sf(data = sgp_boundary, aes()) +
geom_sf(data = res_shp, fill = NA, size = 0.25, color = '#FF8888') +
geom_raster(aes(x, y, fill = value)) +
scale_fill_viridis_c() +
facet_wrap(~variable, ncol = 2) # creating the map
- I think there are strengths and weaknesses to the code used. The use
of 3 common indices, temperature, transport and greenness worked well to
show live-ability as these are some key components to a persons
livelihood, i.e a liveable area cannot be too hot, have no greenery and
must have transport for work/ food etc however, these are not the only
indices that are at play when deciding live-ability as some areas can be
too expensive or lacking in for example shops. Weaknesses lie in the
lack of general public transportation mapping as the MRT is not the only
form of transport in Singapore. Roads could have also been used for
this. Trying to put live-ability to one map is not fully realistic as
every person has their own preferences for live-ability
