Geospatial_PS1
Angad
2025-01-28
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(sf)## Warning: package 'sf' was built under R version 4.4.2## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUElibrary(ggplot2)
library(terra)## Warning: package 'terra' was built under R version 4.4.2## terra 1.8.15# Set working directory
setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 1\\Q1")
# Load Global Dam Watch (GDW) dataset containing dam locations
represas <- st_read("D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 1/Q1/GDW_reservoirs_v1_0.shp")## Reading layer `GDW_reservoirs_v1_0' from data source
## `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 1\Q1\GDW_reservoirs_v1_0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 35295 features and 71 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -159.4638 ymin: -45.88067 xmax: 178.1558 ymax: 70.39624
## Geodetic CRS: WGS 84# Filter dams located in South Africa and constructed before the end of 2002
represas_sa <- represas[represas$COUNTRY == "South Africa" & represas$YEAR_DAM <= 2002, ]
# Load GADM dataset containing district boundaries for South Africa
distritos <- st_read("gadm41_ZAF_3.shp")## Reading layer `gadm41_ZAF_3' from data source
## `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 1\Q1\gadm41_ZAF_3.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 234 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 16.45189 ymin: -34.83514 xmax: 32.89125 ymax: -22.12503
## Geodetic CRS: WGS 84# Load Digital Elevation Model (DEM) - Terrain elevation data
srtm <- rast("South_Africa_SRMT30meters.tif")
# Reduce the resolution of the DEM raster
srtm_reducido <- aggregate(srtm, fact = 5, fun = mean)## |---------|---------|---------|---------|========================================= # Calculate terrain slope in degrees from the DEM
slope <- terrain(srtm_reducido, v = "slope", unit = "degrees")
# Extract the average slope value for each district
gradiente_por_distrito <- extract(slope, vect(distritos), fun = mean, na.rm = TRUE)
# Assign the calculated average slope to the district shapefile
distritos$avg_slope <- gradiente_por_distrito[, 2]
# plot dam locations and terrain slope by district
ggplot() +
# districts, colored by average slope
geom_sf(data = distritos, aes(fill = avg_slope), color = "black") +
# dam locations
geom_sf(data = represas_sa, aes(color = "Dam location"), size = 1) +
# Inverted color scale for slope
scale_fill_viridis_c(option = "plasma", name = "Average Slope", direction = -1) +
# Color scale for dams
scale_color_manual(values = c("Dam location" = "blue"), name = "Legend") +
# Labels and title
labs(title = "Dam Locations and Terrain Slope by District",
caption = "Source: GADM, RCMRD & GDW") 
#Set working directory
# setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 3")
# Load necessary packages
library(sf)
library(geobr)## Warning: package 'geobr' was built under R version 4.4.2library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(readxl)
library(tidyr)##
## Attaching package: 'tidyr'## The following object is masked from 'package:terra':
##
## extractlibrary(lwgeom)## Warning: package 'lwgeom' was built under R version 4.4.2## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.12.2, PROJ 9.4.1##
## Attaching package: 'lwgeom'## The following object is masked from 'package:sf':
##
## st_perimeterlibrary(sfheaders)## Warning: package 'sfheaders' was built under R version 4.4.2# Question 3
# Read the .xls file with the Brazilian population for 2009
df <- read_excel("D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 3/UF_Municipio.xls", skip = 4, col_names = TRUE)## New names:
## • `COD` -> `COD...2`
## • `COD` -> `COD...3`Including Plots
You can also embed plots, for example:
# Get the geographical data from the package geobr
state <- read_state(year = 2010)## Using year/date 2010meso <- read_meso_region(year = 2010)## Using year/date 2010region <- read_region(year = 2010)## Using year/date 2010# Convert ESTIMADA (=population) into a numerical variable and handle missing values
# by setting them to 0
df <- df %>%
mutate(ESTIMADA = as.numeric(ESTIMADA)) %>%
mutate(ESTIMADA = replace_na(ESTIMADA, 0))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `ESTIMADA = as.numeric(ESTIMADA)`.
## Caused by warning:
## ! NAs introduced by coercion# Group by the country code
df_grouped <- df %>%
group_by(COD...2) %>%
summarise(
Total_Estimada = sum(ESTIMADA, na.rm = TRUE),
Count = n()
)
# Remove the last row of df_merged - in df the last row was the total
# and it would have been considered another meso region
df_grouped <- df_grouped %>%
slice(-n())
# Merge the geographical data with the population data
df_merged <- state %>%
left_join(df_grouped, by = c("code_state" = "COD...2"))
# Create the column pop_share
df_merged <- df_merged %>%
mutate(pop_share = Total_Estimada / sum(Total_Estimada) * 100)
# Aggregate regions as the author's did to get the correct boundaries
region <- region %>%
mutate(
name_region = case_when(
name_region %in% c("Centro Oeste", "Norte") ~ 2,
TRUE ~ 1
)
) %>%
dplyr::select(name_region, geom) %>%
group_by(name_region) %>%
summarise(geom = st_union(geom)) %>%
sfheaders::sf_remove_holes()
# Final plot
ggplot(data = df_merged) +
geom_sf(aes(fill = pop_share), color = "white") +
scale_fill_gradientn(
colors = c("lightblue", "blue", "darkblue"),
name = "Pop. Share",
limits = c(0, 25)
) +
geom_sf(data = region, aes(color = factor(name_region)), fill = NA, linewidth = 1) +
scale_color_manual(values = c("black", "red")) +
guides(color = "none")
# Set working directory
setwd("D:\\BSE Semester 2\\Geospatial Data Science and Economic Spatial Models\\Assignments\\Assign 1\\Question 4")# Question 4
# Step 1 Import libraries
library(sf) # simple features' library
library(spData) # library of spatial datasets## Warning: package 'spData' was built under R version 4.4.2## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(tidyverse) # dplyr, ggplot## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ lubridate 1.9.3 ✔ stringr 1.5.1
## ✔ purrr 1.0.2 ✔ tibble 3.2.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors# Step 2 - Load shape files. Note have used different data source to create the map
# Vietnam map, data source - Link : https://gadm.org/download_country.html#google_vignette
sf.vietnam <- st_read('D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 4/gadm41_VNM_shp/gadm41_VNM_0.shp')## Reading layer `gadm41_VNM_0' from data source
## `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 4\gadm41_VNM_shp\gadm41_VNM_0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 102.1446 ymin: 8.381355 xmax: 109.4692 ymax: 23.39269
## Geodetic CRS: WGS 84# Vietnam Roads, data source - Link : https://data.humdata.org/dataset/viet-nam-roads
# Note data is available for 2015
sf.roads <- st_read('D:/BSE Semester 2/Geospatial Data Science and Economic Spatial Models/Assignments/Assign 1/Question 4/vnm_rdsl_2015_osm/vnm_rdsl_2015_OSM.shp')## Reading layer `vnm_rdsl_2015_OSM' from data source
## `D:\BSE Semester 2\Geospatial Data Science and Economic Spatial Models\Assignments\Assign 1\Question 4\vnm_rdsl_2015_osm\vnm_rdsl_2015_OSM.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 118138 features and 5 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 102.1562 ymin: 8.596338 xmax: 109.4536 ymax: 23.37756
## Geodetic CRS: WGS 84# Geo map initial
# ggplot() + geom_sf(data = sf.vietnam) + geom_sf(data= sf.roads, aes(color=type))
# unique(sf.roads$type)
# Step 3 - Club and reduce road types
sf.roads <- sf.roads %>%
mutate(
road_group = case_when(
type %in% c("motorway", "motorway_link", "") ~ "Freeway",
type %in% c("primary", "secondary", "tertiary") ~ "Major Roads",
type %in% c("residential", "service", "track") ~ "Minor Roads",
TRUE ~ "Other"
)
)
# Define custom colors for each road group
road_colors <- c(
"Highways" = "red",
"Major Roads" = "blue",
"Minor Roads" = "green",
"Other" = "gray"
)
# Step 4 - Geo map with custom road groups and colors
ggplot() +
geom_sf(data = sf.vietnam, fill = "lightyellow", color = "black") + # Vietnam map
geom_sf(data = sf.roads, aes(color = road_group)) + # Roads with grouped types
scale_color_manual(values = road_colors, name = "Road Types") + # Custom colors
theme_minimal() +
labs(title = "Road Map of Vietnam (2015)", subtitle = "Road Types Categorized by Function") +
theme(legend.position = "bottom", legend.title = element_text(size = 12))