Loading The data
library(sf)
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
library(tmap)
library(pals)
library(cartography)
## This project is in maintenance mode.
## Core functionalities of `cartography` can be found in `mapsf`.
## https://riatelab.github.io/mapsf/
library(aspace)
## Loading required package: splancs
## Loading required package: sp
##
## Spatial Point Pattern Analysis Code in S-Plus
##
## Version 2 - Spatial and Space-Time analysis
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following object is masked from 'package:splancs':
##
## zoom
## The following objects are masked from 'package:base':
##
## format.pval, units
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::src() masks Hmisc::src()
## ✖ dplyr::summarize() masks Hmisc::summarize()
## ✖ dplyr::tribble() masks tidyr::tribble(), tibble::tribble(), splancs::tribble()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.univar
## spatstat.univar 3.1-2
## Loading required package: spatstat.geom
## spatstat.geom 3.3-6
## Loading required package: spatstat.random
## spatstat.random 3.3-3
## Loading required package: spatstat.explore
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## spatstat.explore 3.4-2
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.3-5
## Loading required package: spatstat.linnet
## spatstat.linnet 3.2-5
##
## spatstat 3.3-2
## For an introduction to spatstat, type 'beginner'
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
School <- st_read("/Users/dhkdtjdud/Desktop/School/上課/113-2/[必修]空間分析/空分期末報告/SCHOOL-3/SCHOOL.shp",options = "encoding=Big5")
## options: encoding=Big5
## Reading layer `SCHOOL' from data source
## `/Users/dhkdtjdud/Desktop/School/上課/113-2/[必修]空間分析/空分期末報告/SCHOOL-3/SCHOOL.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 148 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 297078.6 ymin: 2763290 xmax: 312516.7 ymax: 2784542
## Projected CRS: TWD97 / TM2 zone 121
Tp_vill <- st_read("/Users/dhkdtjdud/Desktop/School/上課/113-2/[必修]空間分析/空分期末報告/Taipei_Vill/Taipei_Vill.shp", options = "encoding=Big5")
## options: encoding=Big5
## Reading layer `Taipei_Vill' from data source
## `/Users/dhkdtjdud/Desktop/School/上課/113-2/[必修]空間分析/空分期末報告/Taipei_Vill/Taipei_Vill.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 456 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 296094.4 ymin: 2761518 xmax: 317198.9 ymax: 2789180
## Projected CRS: TWD97 / TM2 zone 121
Tp_ff <- st_read("/Users/dhkdtjdud/Desktop/School/上課/113-2/[必修]空間分析/空分期末報告/Tpe_Fastfood/Tpe_Fastfood.shp", options = "encoding=Big5")
## options: encoding=Big5
## Reading layer `Tpe_Fastfood' from data source
## `/Users/dhkdtjdud/Desktop/School/上課/113-2/[必修]空間分析/空分期末報告/Tpe_Fastfood/Tpe_Fastfood.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 297198.9 ymin: 2763885 xmax: 312205.7 ymax: 2781148
## Projected CRS: TWD97 / TM2 zone 121
School <- st_transform(School, 3826)
Tp_vill <- st_transform(Tp_vill, 3826)
Tp_ff <- st_transform(Tp_ff, 3826)
study_area <- st_union(Tp_vill)
grid <- st_make_grid(Tp_vill, cellsize = 500, square = TRUE)
grid_sf <- st_sf(geometry = grid)
grid_sf$grid_id <- 1:nrow(grid_sf)
ff_joined <- st_join(Tp_ff, grid_sf, join = st_within, left = FALSE)
ff_count <- ff_joined %>%
st_drop_geometry() %>%
group_by(grid_id) %>%
summarise(store_count = n())
grid_sf <- left_join(grid_sf, ff_count, by = "grid_id")
grid_sf$store_count[is.na(grid_sf$store_count)] <- 0
school_joined <- st_join(School, grid_sf, join = st_within, left = FALSE)
school_count <- school_joined %>%
st_drop_geometry() %>%
group_by(grid_id) %>%
summarise(school_count = n())
grid_sf <- left_join(grid_sf, school_count, by = "grid_id")
grid_sf$school_count[is.na(grid_sf$school_count)] <- 0
nb <- poly2nb(grid_sf, queen = FALSE)
lw <- nb2listw(nb, style = "B", zero.policy = TRUE)
grid_sf$store_z <- scale(grid_sf$store_count)
gi_store <- localG(as.vector(grid_sf$store_z), lw, zero.policy = TRUE)
grid_sf$Gi_store <- as.numeric(gi_store)
grid_sf$Gi_store_sig <- ifelse(grid_sf$Gi_store > 1.96, "Hot Spot",
ifelse(grid_sf$Gi_store < -1.96, "Cold Spot", "Not Significant"))
grid_sf$school_z <- scale(grid_sf$school_count)
gi_school <- localG(as.vector(grid_sf$school_z), lw, zero.policy = TRUE)
grid_sf$Gi_school <- as.numeric(gi_school)
grid_sf$Gi_school_sig <- ifelse(grid_sf$Gi_school > 1.96, "Hot Spot",
ifelse(grid_sf$Gi_school < -1.96, "Cold Spot", "Not Significant"))
tmap_mode("plot")
## ℹ tmap mode set to "plot".
# 地圖 1：速食店熱區
map1 <- tm_shape(grid_sf) +
tm_fill("Gi_store_sig", palette = c("red", "blue", "gray"),
title = "Fast Food\nGi* Hotspot",
labels = c("Hot Spot (p<0.05)", "Cold Spot (p<0.05)", "Not Significant")) +
tm_borders(col = "black", lwd = 0.1) +
tm_layout(title = "Getis-Ord Gi* Hot Spot Analysis(Fast Food)", legend.outside = TRUE)
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values'), 'labels' to
## fill.scale = tm_scale(<HERE>).[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
# 地圖 2：學校熱區
map2 <- tm_shape(grid_sf) +
tm_fill("Gi_school_sig", palette = c("red", "blue", "gray"),
title = "School\nGi* Hotspot",
labels = c("Hot Spot (p<0.05)", "Cold Spot (p<0.05)", "Not Significant")) +
tm_borders(col = "black", lwd = 0.1) +
tm_layout(title = "Getis-Ord Gi* Hot Spot Analysis(School)", legend.outside = TRUE)
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
# 顯示地圖
tmap_arrange(map1, map2)
## Warning: labels do not have the same length as levels, so they are repeated
## Warning: labels do not have the same length as levels, so they are repeated
## Warning: labels do not have the same length as levels, so they are repeated
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
## Warning: labels do not have the same length as levels, so they are repeated
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
