Homework 1
Student number:A0326906Y Name:Wang Maoye
2025-09-29
Question 1
a.
library(spatstat.data)
library(spatstat.geom)
library(spatstat.explore)
library(spatstat)
data(longleaf)
plot(longleaf, main = "Raw Spatial Point Pattern of Longleaf Trees")
- Monte Carlo test based on G function
G_test <- envelope(longleaf,
fun = Gest, # G function
nsim = 199, # number of simulations
global = FALSE) # equal-tailed envelope## Generating 199 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34
## .36.38.40.42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74
## .76.78.80.82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114
## .116.118.120.122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154
## .156.158.160.162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194
## .196.198
## 199.
##
## Done.plot(G_test, main = "G Function Monte Carlo Test for CSR")
For the G function, most of its values lie above the 95% confidence band, suggesting a clustering tendency in the data.
- Monte Carlo test based on F function
F_test <- envelope(longleaf,
fun = Fest, # F function
nsim = 199,
global = FALSE) # equal-tailed envelope## Generating 199 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34
## .36.38.40.42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74
## .76.78.80.82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114
## .116.118.120.122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154
## .156.158.160.162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194
## .196.198
## 199.
##
## Done.plot(F_test, main = "F Function Monte Carlo Test for CSR")
In the case of the F function, the majority of its values fall below the 95% confidence band, which similarly points to clustering in the data.
- Quadrat Test (4x4 Grid)
Q_counts <- quadratcount(longleaf, nx = 4, ny = 4) # Count points in each subregion
Q_chisq <- chisq.test(as.vector(Q_counts)) # Chi-square test
print(Q_chisq) # Display test result in PDF##
## Chi-squared test for given probabilities
##
## data: as.vector(Q_counts)
## X-squared = 103.78, df = 15, p-value = 2.495e-15plot(longleaf, main = "Quadrat Test (4x4 Grid)")
plot(Q_counts, add = TRUE, col = "red") # Display counts on the grid
Q_chisq <- chisq.test(as.vector(Q_counts)) # Chi-square test
print(Q_chisq) # Display test result in PDF##
## Chi-squared test for given probabilities
##
## data: as.vector(Q_counts)
## X-squared = 103.78, df = 15, p-value = 2.495e-15The quadrat test with a 4 × 4 subdivision yields a p-value of 4.99 × 10−15, far below the 0.05 significance threshold. Hence, the null hypothesis of CSR is rejected. When visualizing the longleaf dataset, the spatial points appear clustered, consistent with the statistical test outcome.
Conclusion: The three test results agree with the spatial point pattern observed from running plot(longleaf).
b.
Plot spatial point patterns for each group
# Split longleaf into adult (dbh >= 30cm) and juvenile (dbh < 30cm) trees
adult <- longleaf[marks(longleaf) >= 30] # Adult trees
juvenile <- longleaf[marks(longleaf) < 30] # Juvenile trees
# Plot spatial point patterns for each group
plot(adult, main = "Spatial Point Pattern of Adult Trees (dbh >= 30cm)") 
plot(juvenile, main = "Spatial Point Pattern of Juvenile Trees (dbh < 30cm)") 
The spatial point pattern plots indicate that adult trees are relatively dispersed, whereas juvenile trees exhibit a clear tendency to cluster.
- G Function Monte Carlo Test for each group
# Adult trees
G_adult <- envelope(adult, fun = Gest, nsim = 199, global = FALSE) # Compute G function envelope## Generating 199 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34
## .36.38.40.42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74
## .76.78.80.82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114
## .116.118.120.122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154
## .156.158.160.162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194
## .196.198
## 199.
##
## Done.plot(G_adult, main = "G Function Envelope - Adult Trees")
# Juvenile trees
G_juvenile <- envelope(juvenile, fun = Gest, nsim = 199, global = FALSE)## Generating 199 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34
## .36.38.40.42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74
## .76.78.80.82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114
## .116.118.120.122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154
## .156.158.160.162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194
## .196.198
## 199.
##
## Done.plot(G_juvenile, main = "G Function Envelope - Juvenile Trees")
The G function results show that nearly all values for juvenile trees lie above the 95% confidence band, indicating strong clustering, while adult trees display values closer to the confidence band, suggesting less pronounced clustering.
- F Function Monte Carlo Test for each group
# Adult trees
F_adult <- envelope(adult, fun = Fest, nsim = 199, global = FALSE) # Compute F function envelope## Generating 199 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34
## .36.38.40.42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74
## .76.78.80.82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114
## .116.118.120.122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154
## .156.158.160.162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194
## .196.198
## 199.
##
## Done.plot(F_adult, main = "F Function Envelope - Adult Trees")
# Juvenile trees
F_juvenile <- envelope(juvenile, fun = Fest, nsim = 199, global = FALSE)## Generating 199 simulations of CSR ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34
## .36.38.40.42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74
## .76.78.80.82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114
## .116.118.120.122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154
## .156.158.160.162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194
## .196.198
## 199.
##
## Done.plot(F_juvenile, main = "F Function Envelope - Juvenile Trees")
The F function results reveal that almost all points for juvenile trees fall below the 95% confidence band, confirming significant clustering, whereas adult trees show minimal deviation from the band, indicating weaker clustering.
- Quadrat Test for each group
# Adult Trees
# -----------------------------
# Count points in each 4x4 quadrat
Q_counts_adult <- quadratcount(adult, nx = 4, ny = 4)
# Perform chi-squared test against CSR expectation
Q_chisq_adult <- chisq.test(as.vector(Q_counts_adult))
print(Q_chisq_adult)##
## Chi-squared test for given probabilities
##
## data: as.vector(Q_counts_adult)
## X-squared = 63.465, df = 15, p-value = 6.337e-08# Plot adult trees and overlay quadrat counts in red
plot(adult, main = "Quadrat Test - Adult Trees (4x4 Grid)")
plot(Q_counts_adult, add = TRUE, col = "red") # Overlay counts on the grid
# Juvenile Trees
# -----------------------------
# Count points in each 4x4 quadrat
Q_counts_juvenile <- quadratcount(juvenile, nx = 4, ny = 4)
# Perform chi-squared test against CSR expectation
Q_chisq_juvenile <- chisq.test(as.vector(Q_counts_juvenile))
print(Q_chisq_juvenile) ##
## Chi-squared test for given probabilities
##
## data: as.vector(Q_counts_juvenile)
## X-squared = 233.4, df = 15, p-value < 2.2e-16# Plot juvenile trees and overlay quadrat counts in red
plot(juvenile, main = "Quadrat Test - Juvenile Trees (4x4 Grid)")
plot(Q_counts_juvenile, add = TRUE, col = "red") # Overlay counts on the grid
The quadrat test yields a p-value < 2.2 × 10⁻¹⁶ and X² = 233.4 for juvenile trees, compared to p = 6.337 × 10⁻⁸ and X² = 63.465 for adult trees, again highlighting stronger clustering among juvenile trees.
Question 2
a.
library(sf) # For handling spatial objects (sf)
library(ggplot2) # For plotting mapsStep 1: Read the RDS file into R
# Replace the path with your actual file location
file_path <- "C:/Users/dydsw/Desktop/sg_subzones.rds"
# Read the RDS file; it should contain an sf object with 310 subzones
sg_subzones <- readRDS(file_path)
# Verify that the object is indeed an sf object
class(sg_subzones) # Expected: "sf" "data.frame"## [1] "sf" "data.frame"Step 2: Transform the coordinate reference system to EPSG:4326
# Transform to WGS84 geographic coordinates (longitude/latitude)
sg_subzones1 <- st_transform(sg_subzones, crs = 4326)
# Check the new CRS to confirm
st_crs(sg_subzones1)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## 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]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]Step 3: Merge all subzones into a single polygon
# Use st_union to combine all subzones, removing internal boundaries
# 'by_feature = FALSE' ensures we get a single geometry
sg_union <- st_union(sg_subzones1, by_feature = FALSE)
# Convert the unioned geometry back to an sf object for ggplot
sg_union_sf <- st_sf(geometry = sg_union)Step 4: Plot the Singapore map without internal subzone boundaries
# geom_sf() automatically handles sf objects
# fill = "lightblue" sets the fill color
# color = "black" sets the polygon outline
ggplot(sg_union_sf) +
geom_sf(fill = "lightblue", color = "black") +
ggtitle("Singapore Map (Merged Subzones, EPSG:4326)") +
theme_minimal()
b.
# Filter out rows where the 'total' population is NA
sg_subzones2 <- sg_subzones1[!is.na(sg_subzones1$total), ]
# Verify the filtering
summary(sg_subzones2$total)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 3730 11470 15792 137630# Compute log10(total + 1) to avoid log(0)
sg_subzones2$log_total <- log10(sg_subzones2$total + 1)
# Check the first few values
head(sg_subzones2$log_total)## [1] 4.210345 3.645521 0.000000 3.380392 4.082821 4.199508# Use ggplot2 to visualize subzones colored by log-transformed population
ggplot(sg_subzones2) +
geom_sf(aes(fill = log_total), color = "black") + # fill color based on log_total
scale_fill_viridis_c(option = "plasma") + # use perceptually uniform color scale
ggtitle("Singapore Subzones Colored by log10(Total Population + 1)") +
theme_minimal()
c.
# Set the path to the downloaded RDS file
file_path_2015 <- "C:/Users/dydsw/Desktop/sg_subzones_2015.rds"
# Read the RDS file into an sf object
sg.subzones.2015 <- readRDS(file_path_2015)
class(sg.subzones.2015) # Verify that it's an sf object## [1] "sf" "data.frame"# Compute log10(POP + 1) to avoid log(0)
sg.subzones.2015$log_POP <- log10(sg.subzones.2015$POP + 1)
# Use ggplot2 to visualize all subzones, coloring by log_POP
ggplot(sg.subzones.2015) +
geom_sf(aes(fill = log_POP)) + # Fill color mapped to log_POP
scale_fill_viridis_c() + # Continuous color scale
ggtitle("Singapore Subzones 2015 - log_POP") + # Add a title
theme_minimal() # Minimal theme for clarity
d.
# Population density is defined as POP divided by SHAPE_Area
sg.subzones.2015$density <- sg.subzones.2015$POP / sg.subzones.2015$SHAPE_Area
# Order by POP in descending order
top5_pop <- sg.subzones.2015[order(-sg.subzones.2015$POP), ][1:5, ]
# Print the names of the top 5 most populated subzones
cat("Top 5 most populated subzones:\n")## Top 5 most populated subzones:print(top5_pop$SUBZONE_N)## [1] "TAMPINES EAST" "WOODLANDS EAST" "BEDOK NORTH" "TAMPINES WEST"
## [5] "YUNNAN"The top 5 most densely populated subzones are TAMPINES EAST、WOODLANDS EAST、BEDOK NORTH、TAMPINES WEST、YUNNAN
Question 3
a.
# Load required packages
library(sf)
# Read station shapefile as sf object
stations <- st_read("C:/Users/dydsw/Desktop/met_stations/stations.shp")## Reading layer `stations' from data source
## `C:\Users\dydsw\Desktop\met_stations\stations.shp' using driver `ESRI Shapefile'
## Simple feature collection with 63 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 103.6184 ymin: 1.189 xmax: 103.9826 ymax: 1.4582
## Geodetic CRS: WGS 84# Read rainfall CSV file
all_data <- read.csv("C:/Users/dydsw/Desktop/rainfall_201601.csv",
stringsAsFactors = FALSE)
# Check data structures
head(stations) # show first rows of spatial data## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 103.7042 ymin: 1.3087 xmax: 103.8492 ymax: 1.4582
## Geodetic CRS: WGS 84
## code name geometry
## 1 S104 Admiralty POINT (103.7854 1.4439)
## 2 S105 Admiralty (West) POINT (103.7953 1.4582)
## 3 S109 Ang Mo Kio POINT (103.8492 1.3764)
## 4 S86 Boon Lay (East) POINT (103.7205 1.3302)
## 5 S63 Boon Lay (West) POINT (103.7042 1.3275)
## 6 S120 Botanic Garden POINT (103.818 1.3087)head(all_data) # show first rows of rainfall data## stn_name year month day daily_rainfall stn_code
## 1 Yishun 2016 1 1 6.4 S91
## 2 Yishun 2016 1 2 23.8 S91
## 3 Yishun 2016 1 3 1.2 S91
## 4 Yishun 2016 1 4 0.2 S91
## 5 Yishun 2016 1 5 0.0 S91
## 6 Yishun 2016 1 6 0.0 S91b.
library(sf)
met.stations <- st_read("C:/Users/dydsw/Desktop/met_stations/stations.shp", quiet = TRUE)
all_data <- read.csv("C:/Users/dydsw/Desktop/rainfall_201601.csv", stringsAsFactors = FALSE)
# Using aggregate() as suggested in hint
avg_rainfall <- aggregate(daily_rainfall ~ stn_code, data = all_data, FUN = mean, na.rm = TRUE)
colnames(avg_rainfall) <- c("code", "avg_rainfall") # Rename columns for mergingrm = TRUE))
# Perform a left join to keep all stations and add mean_rainfall
met.stations <- merge(met.stations, avg_rainfall, by = "code", all.x = TRUE)
# Check result
head(met.stations) # Display first few rows with average rainfall## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 103.7134 ymin: 1.189 xmax: 103.9007 ymax: 1.4439
## Geodetic CRS: WGS 84
## code name avg_rainfall geometry
## 1 S06 Paya Lebar 2.774194 POINT (103.9007 1.3524)
## 2 S07 Macritchie\nReservoir 5.587097 POINT (103.8338 1.3417)
## 3 S08 Lower Peirce\nReservoir 6.090323 POINT (103.8271 1.3701)
## 4 S101 Jurong (North) 10.155556 POINT (103.7134 1.3505)
## 5 S102 Semakau Island 6.993548 POINT (103.768 1.189)
## 6 S104 Admiralty 5.625806 POINT (103.7854 1.4439)c.
# Load required libraries
# ---------------------------------------------
library(sf)
library(ggplot2)
# Load the shapefile as an sf object
sg_boundary <- st_read("C:/Users/dydsw/Desktop/Singapore_Boundary/Singapore_Boundary.shp", quiet = TRUE)
# Load Singapore boundary shapefile
singapore_boundary <- st_read("C:/Users/dydsw/Desktop/Singapore_Boundary/Singapore_Boundary.shp", quiet = TRUE)
# Plot Singapore boundary with average rainfall at stations
ggplot() +
geom_sf(data = singapore_boundary, fill = "lightgrey", color = "black") + # Singapore map
geom_sf(data = met.stations, aes(color = avg_rainfall), size = 2) + # rainfall points
scale_color_viridis_c(option = "viridis", na.value = "grey50") + # color scale
ggtitle("Average Rainfall of 63 Stations in Jan 2016") +
theme_minimal()
Question 4
a.
library(geodata)
library(terra)
library(sf)
library(ggplot2)
# Download max temperature (tmax) for Singapore, specify path to save data
tmax_sg <- worldclim_country(country = "SGP", var = "tmax", path = "C:/Users/dydsw/Desktop/WorldClimData")
# Average over the 12 monthly layers and save as sg.tmax
sg.tmax <- app(tmax_sg, fun = mean)
# Read Singapore boundary
sg_border <- st_read("C:/Users/dydsw/Desktop/Singapore_Boundary/Singapore_Boundary.shp", quiet = TRUE)
# Crop and mask temperature raster
sg.tmax <- crop(sg.tmax, sg_border)
sg.tmax <- mask(sg.tmax, sg_border)
# Plot
plot(sg.tmax, main = "Average Maximum Temperature in Singapore (°C)")
plot(st_geometry(sg_border), add = TRUE, border = "black")
b.
# Download precipitation (prec) for Singapore
prec_sg <- worldclim_country(country = "SGP", var = "prec", path = "C:/Users/dydsw/Desktop/WorldClimData")
# Average over 12 months and save as sg.prec
sg.prec <- app(prec_sg, fun = mean)
# Crop and mask precipitation raster
sg.prec <- crop(sg.prec, sg_border)
sg.prec <- mask(sg.prec, sg_border)
# Plot
plot(sg.prec, main = "Average Precipitation in Singapore")
plot(st_geometry(sg_border), add = TRUE, border = "black")
c.
# Align rasters: crop to common extent
common_ext <- intersect(ext(sg.tmax), ext(sg.prec))
tmax_crop <- crop(sg.tmax, common_ext)
prec_crop <- crop(sg.prec, common_ext)
# Extract values
tmax_vals <- values(tmax_crop)
prec_vals <- values(prec_crop)
# Combine into a matrix and remove rows with NA
df <- na.omit(cbind(tmax_vals, prec_vals))
# Compute sample correlation
cor(df)## mean mean
## mean 1.0000000 0.3981042
## mean 0.3981042 1.0000000d.
# Read Master Plan subzones shapefile
master_subzones <- st_read("C:/Users/dydsw/Desktop/Singapore_Subzones/MasterPlan2019SubzoneBoundary.shp", quiet = TRUE)
# Extract average precipitation per subzone
avg_prec_values <- terra::extract(sg.prec, vect(master_subzones), fun = mean, na.rm = TRUE, weights = FALSE)
wtd_prec_values <- terra::extract(sg.prec, vect(master_subzones), fun = mean, na.rm = TRUE, weights = TRUE)
# Add extracted values as new columns
master_subzones$avg_prec <- avg_prec_values[,2] # column 2 contains the extracted values
master_subzones$wtd_avg_prec <- wtd_prec_values[,2]
# Plot average precipitation
p1 <- ggplot(master_subzones) +
geom_sf(aes(fill = avg_prec)) +
scale_fill_viridis_c() +
ggtitle("Average Precipitation per Subzone") +
theme_bw()
# Plot weighted average precipitation
p2 <- ggplot(master_subzones) +
geom_sf(aes(fill = wtd_avg_prec)) +
scale_fill_viridis_c() +
ggtitle("Weighted Average Precipitation per Subzone") +
theme_bw()
# Arrange plots side by side
library(patchwork)
p1 + p2 + plot_layout(widths = c(1, 1))