final

rm(list = ls())
library(spdep)

## 載入需要的套件：spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## 載入需要的套件：sf

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

library(sf)
library(tmap)
library(dplyr)

## 
## 載入套件：'dplyr'

## 下列物件被遮斷自 'package:stats':
## 
##     filter, lag

## 下列物件被遮斷自 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

setwd("C://Users/user/Desktop/spacial_anl/")
school=st_read("SCHOOL.shp")

## Reading layer `SCHOOL' from data source 
##   `C:\Users\user\Desktop\spacial_anl\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("Taipei_Vill.shp",options="ENCODING=BIG5")

## options:        ENCODING=BIG5 
## Reading layer `Taipei_Vill' from data source 
##   `C:\Users\user\Desktop\spacial_anl\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_food=st_read("Tpe_Fastfood.shp")

## Reading layer `Tpe_Fastfood' from data source 
##   `C:\Users\user\Desktop\spacial_anl\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, crs = 3826)
tp_vill<- st_transform(tp_vill, crs = 3826)
tp_food <- st_transform(tp_food, crs = 3826)

500公尺網格

grid <- st_make_grid(
  tp_vill,
  cellsize = 500,
  what = "polygons",
  square = TRUE
) %>%
  st_sf() %>%
  mutate(grid_id = row_number())

統計每個格子內的速食店數

# 交集統計-速食店數
grid_count <- st_join(grid, tp_food, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(faststore_count = n())

grid_count$faststore_count[is.na(grid_count$faststore_count)] <- 0

qeen

grid.nb <- poly2nb(grid_count, queen = TRUE)  #角落接觸也算
grid.nb.w <- nb2listw(grid.nb, style = "W", zero.policy = T)

grid.nb.in=include.self(grid.nb)

grid.nb.w.in=nb2listw(grid.nb.in)

faststore_count=grid_count$faststore_count

Gi=localG(faststore_count, grid.nb.w.in)

Bonferroni多重檢定修正

n=nrow(grid)

grid$Gi=localG(faststore_count,grid.nb.w.in)
grid$colG="grey99"

#Bonferroni多重檢定修正
grid$colG[grid$Gi>=qnorm(1-0.05/n)]="red" 

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(c())

district_colors <- c(
  "#1b9e77", "#d95f02", "#7570b3", "#e7298a",
  "#66a61e", "#e6ab02", "#a6761d", "#666666",
  "#a6cee3", "#fb9a99", "#fdbf6f", "#b2df8a"
)
tm_shape(grid) +
  tm_fill("colG", title = "速食店熱區", palette = c("grey99", "green")) +
  tm_borders(col = "lightblue") +
  
  tm_shape(tp_vill)+tm_borders("grey")+

  tm_shape(tp_district) +
  tm_fill("TOWN", alpha = 0.3, title = "行政區", palette = district_colors) +   # 加透明度，不蓋住底圖
  tm_borders(lwd = 1, col = "black") +

  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("速食店500m網格熱區")

## 

## ── 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') 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_fill()`: use `fill_alpha` instead of `alpha`.
## [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>)'

由QGIS軟體，批次轉換台北市高中職、國中小的地址為有空間屬性的sf檔案

開啟台北市高中職、國中小資料

tp_school=st_read("C:/Users/user/Downloads/tp_school.shp",options="ENCODING=BIG5")

## options:        ENCODING=BIG5 
## Reading layer `tp_school' from data source `C:\Users\user\Downloads\tp_school.shp' using driver `ESRI Shapefile'
## Simple feature collection with 282 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 121.4671 ymin: 24.97251 xmax: 121.624 ymax: 25.16832
## Geodetic CRS:  WGS 84

tp_school <- st_transform(tp_school, crs = 3826)

tp_school_lyr=tm_shape(tp_school)+tm_dots("#00008B")

tp_lyr=tm_shape(tp_vill)+tm_polygons("grey90")+tm_borders("grey")+tm_shape(tp_district)+tm_borders("black",lwd=1.5)
tp_lyr+tm_layout(frame=F)+tp_school_lyr

發現有一筆資料在台北市外，因此要找到它並修正座標。

找到台北市外的點，並更正座標

intersections <- st_intersects(tp_school, tp_vill)

# 找出沒與任何 polygon 相交的點（也就是在外面的）
outside_idx <- which(lengths(intersections) == 0)

outside_points <- tp_school[outside_idx, ]
outside_points[,c("schoolname","address")]

#  「仁愛國小 臺北市大安區敦安里安和路1段60號」資料點在台北市外

tp_lyr+tm_shape(outside_points)+tm_dots("red",size=0.5)

tp_school=st_transform(tp_school,crs=4326) #經緯度

# 座標來源：google map
new_geom <- st_sfc(st_point(c(121.55216829140444, 25.03605419284852)), crs = 4326) 

# 更新第 23 筆資料的 geometry 欄位
st_geometry(tp_school)[23] <- new_geom

tp_school=st_transform(tp_school,3826)

tp_school_lyr=tm_shape(tp_school)+tm_dots("#00008B")
tp_lyr+tm_layout(frame=F)+tp_school_lyr+tm_shape(tp_school[23,])+tm_dots("red",size=0.5)

更正座標完畢。

速食店熱區在哪些區和里，包含哪些學校

公園點資料

park <- read.csv("C:/Users/user/Downloads/park.csv")

# 轉換為 sf 物件，指定 TW97 投影 EPSG:3826
park<- st_as_sf(park, coords = c("TW97座標X", "TW97座標Y"), crs = 3826)

# 檢查
print(park)

## Simple feature collection with 517 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 296961.9 ymin: 2763622 xmax: 323503.4 ymax: 2782707
## Projected CRS: TWD97 / TM2 zone 121
## First 10 features:
##    ID 行政區   公園名稱 無障礙出入口名稱 人行道寬度 坡度
## 1   1 萬華區   青年公園        1號出入口       1.80    8
## 2   5 萬華區 和平青草園        1號出入口       2.20    5
## 3   7 萬華區   艋舺公園        1號出入口      11.20    6
## 4   9 萬華區   日善公園        1號出入口       3.45    6
## 5  12 萬華區   壽德公園        1號出入口       2.10    8
## 6  13 萬華區   桂林公園        1號出入口      11.60    4
## 7  14 萬華區   華西公園        1號出入口       4.30    3
## 8  16 萬華區   青山公園        1號出入口       1.80    7
## 9  17 萬華區   西園公園        1號出入口       2.00    6
## 10 20 萬華區   長沙公園        1號出入口       3.20    3
##                    geometry
## 1  POINT (301105.5 2768260)
## 2  POINT (300167.6 2769519)
## 3  POINT (300398.4 2769794)
## 4  POINT (300101.6 2768764)
## 5  POINT (300569.8 2768596)
## 6    POINT (300139 2770199)
## 7  POINT (300305.8 2770145)
## 8  POINT (300295.6 2770249)
## 9  POINT (300232.2 2768864)
## 10 POINT (300465.3 2770502)

tp_lyr+tm_layout(frame=F)+tm_shape(park)+tm_dots("#008b45")

district_park=st_intersection(park,tp_district,join=within)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

district_park=district_park %>% group_by(TOWN) %>% summarise(count=n())

#使資料變df，來做left_join
district_park_df=st_drop_geometry(district_park)
tp_district_df=as.data.frame(tp_district)

district_park=left_join(district_park_df,tp_district_df,by="TOWN")
district_park=st_as_sf(district_park)

##密度
district_park$area=st_area(district_park)
district_park$dens=district_park$count/district_park$area

tm_shape(district_park) +
  tm_fill("dens", palette = "YlGn", title = "每區公園密度(1/km^2)",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("每區公園密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"
## Multiple palettes called "yl_gn" found: "brewer.yl_gn", "matplotlib.yl_gn". The first one, "brewer.yl_gn", is returned.

速食店資料

fastfood_lyr=tm_shape(tp_food)+tm_dots("red")+tm_layout(frame=F)
tp_lyr+fastfood_lyr

district_fast=st_intersection(tp_food,tp_district,join=within)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

district_fast=district_fast %>% group_by(TOWN) %>% summarise(count=n())

#使資料變df，來做left_join
district_fast_df=st_drop_geometry(district_fast)
tp_district_df=as.data.frame(tp_district)

district_fast=left_join(district_fast_df,tp_district_df,by="TOWN")
district_fast=st_as_sf(district_fast)

##密度
district_fast$area=st_area(district_fast)
district_fast$dens=district_fast$count/district_fast$area

tm_shape(district_fast) +
  tm_fill("dens", palette = "OrRd", title = "每區速食店密度(1/km^2)",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("每區速食店密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named
## "brewer.or_rd"
## Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

飲料店資料(原始資料包含台北市、新北市)

drink=st_read("C:/Users/user/Downloads/drinks_clean/drinks_clean.gpkg",options="ENCODING=BIG5")

## options:        ENCODING=BIG5

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 6: driver GPKG does not support open option ENCODING

## Reading layer `drinks_clean' from data source 
##   `C:\Users\user\Downloads\drinks_clean\drinks_clean.gpkg' using driver `GPKG'
## Simple feature collection with 1655 features and 24 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 121.4127 ymin: 24.93438 xmax: 121.7104 ymax: 25.23992
## Geodetic CRS:  WGS 84

drink=st_transform(drink,crs=3826)
#保留在臺北市內的飲料店資料
drink <- st_join(drink, tp_vill, join = st_within, left = FALSE)


drink_lyr=tm_shape(drink)+tm_dots("orange")+tm_layout(frame=F)
tp_lyr+drink_lyr

district_drink=st_intersection(drink,tp_district,join=within)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

district_drink=district_drink %>% group_by(TOWN) %>% summarise(count=n())

#使資料變df，來做left_join
district_drink_df=st_drop_geometry(district_drink)
tp_district_df=as.data.frame(tp_district)

district_drink=left_join(district_drink_df,tp_district_df,by="TOWN")
district_drink=st_as_sf(district_drink)

##密度
district_drink$area=st_area(district_drink)
district_drink$dens=district_drink$count/district_drink$area

tm_shape(district_drink) +
  tm_fill("dens", palette = "OrRd", title = "每區飲料店密度(1/km^2)",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("每區飲料店密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named
## "brewer.or_rd"
## Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

游泳池

tp_pool=st_read("C:/Users/user/Downloads/tp_pools.shp",options="ENCODING=BIG5")

## options:        ENCODING=BIG5 
## Reading layer `tp_pools' from data source `C:\Users\user\Downloads\tp_pools.shp' using driver `ESRI Shapefile'
## Simple feature collection with 78 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 121.4671 ymin: 24.98698 xmax: 121.624 ymax: 25.15131
## Geodetic CRS:  WGS 84

tp_pool=st_transform(tp_pool,crs=3826)

pool_lyr=tm_shape(tp_pool)+tm_dots("#4292c6")+tm_layout(frame=F)
tp_lyr+pool_lyr

district_pool=st_intersection(tp_pool,tp_district,join=within)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

district_pool=district_pool %>% group_by(TOWN) %>% summarise(count=n())

#使資料變df，來做left_join
district_pool_df=st_drop_geometry(district_pool)
tp_district_df=as.data.frame(tp_district)

district_pool=left_join(district_pool_df,tp_district_df,by="TOWN")
district_pool=st_as_sf(district_pool)

##密度
district_pool$area=st_area(district_pool)
district_pool$dens=district_pool$count/district_pool$area

tm_shape(district_pool) +
  tm_fill("dens", palette = "YlGn", title = "每區泳池密度(1/km^2)",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("每區泳池密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"
## Multiple palettes called "yl_gn" found: "brewer.yl_gn", "matplotlib.yl_gn". The first one, "brewer.yl_gn", is returned.

免費運動場地

sport=st_read("C:/Users/user/Downloads/free_open_sport_tp.gpkg")

## Reading layer `free_open_sport_tp' from data source 
##   `C:\Users\user\Downloads\free_open_sport_tp.gpkg' using driver `GPKG'
## Simple feature collection with 409 features and 23 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 121.4632 ymin: 24.981 xmax: 121.6245 ymax: 25.1689
## Geodetic CRS:  WGS 84

sport=st_transform(sport,crs=3826)

sport_lyr=tm_shape(sport)+tm_dots("purple")+tm_layout(frame=F)
tp_lyr+sport_lyr

district_sport=st_intersection(sport,tp_district,join=within)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

district_sport=district_sport %>% group_by(TOWN) %>% summarise(count=n())

#使資料變df，來做left_join
district_sport_df=st_drop_geometry(district_sport)
tp_district_df=as.data.frame(tp_district)

district_sport=left_join(district_sport_df,tp_district_df,by="TOWN")
district_sport=st_as_sf(district_sport)

##密度
district_sport$area=st_area(district_sport)
district_sport$dens=district_sport$count/district_sport$area

tm_shape(district_sport) +
  tm_fill("dens", palette = "YlGn", title = "每區免費運動場密度(1/km^2)",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("每區免費運動場密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGn" is named
## "brewer.yl_gn"
## Multiple palettes called "yl_gn" found: "brewer.yl_gn", "matplotlib.yl_gn". The first one, "brewer.yl_gn", is returned.

交集統計-公園數、游泳池數、免費運動場數、速食店數、飲料店數

park_count <- st_join(grid, park, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(park_count = n())

park_count$park_count[is.na(grid_count$park_count)] <- 0

## Warning: Unknown or uninitialised column: `park_count`.

pool_count <- st_join(grid, tp_pool, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(pool_count = n())

pool_count$pool_count[is.na(grid_count$pool_count)] <- 0

## Warning: Unknown or uninitialised column: `pool_count`.

sport_count <- st_join(grid, sport, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(sport_count = n())

sport_count$sport_count[is.na(grid_count$sport_count)] <- 0

## Warning: Unknown or uninitialised column: `sport_count`.

drink_count <- st_join(grid, drink, join = st_contains) %>%
  group_by(grid_id) %>%
  summarise(drink_count = n())

drink_count$drink_count[is.na(grid_count$drink_count)] <- 0

## Warning: Unknown or uninitialised column: `drink_count`.

#sf轉df

park_count_df=st_drop_geometry(park_count)
pool_count_df=st_drop_geometry(pool_count)
sport_count_df=st_drop_geometry(sport_count)
faststore_count_df=st_drop_geometry(grid_count)
drink_count_df=st_drop_geometry(drink_count)


grid_count <- grid %>%
  left_join(park_count_df, by = "grid_id") %>%
  left_join(pool_count_df, by = "grid_id") %>%
  left_join(sport_count_df,by="grid_id") %>%
  left_join(faststore_count_df, by = "grid_id") %>%
  left_join(drink_count_df,by="grid_id") 

公園網格熱區

park_count=as.numeric(park_count$park_count)
grid.park<-grid
grid.park$Gi_park=localG(park_count,grid.nb.w.in)
grid.park$colG_park="grey99"

#Bonferroni多重檢定修正
grid.park$colG_park[grid.park$Gi_park>=qnorm(1-0.05/n)]="green" 

公園網格熱區與行政區疊圖

district_colors <- c(
  "#1b9e77", "#d95f02", "#7570b3", "#e7298a",
  "#66a61e", "#e6ab02", "#a6761d", "#666666",
  "#a6cee3", "#fb9a99", "#fdbf6f", "#b2df8a"
)
tm_shape(grid.park) +
  tm_fill("colG_park", title = "公園熱區", palette = c("grey99", "green")) +
  tm_borders(col = "lightblue") +
  
  tm_shape(tp_vill)+tm_borders("grey")+

  tm_shape(tp_district) +
  tm_fill("TOWN", alpha = 0.3, title = "行政區", palette = district_colors) +   # 加透明度，不蓋住底圖
  tm_borders(lwd = 1, col = "black") +

  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("公園500m網格熱區")

## 

## ── 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') 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_fill()`: use `fill_alpha` instead of `alpha`.
## [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>)'

游泳池網格熱區

pool_count=as.numeric(pool_count$pool_count)
grid.pool<-grid
grid.pool$Gi_pool=localG(pool_count,grid.nb.w.in)
grid.pool$colG_pool="grey99"

#Bonferroni多重檢定修正
grid.pool$colG_pool[grid.pool$Gi_pool>=qnorm(1-0.05/n)]="skyblue" 

游泳池網格熱區與行政區疊圖

tm_shape(grid.pool) +
  tm_fill("colG_pool", title = "游泳池熱區", palette = c("grey99", "green")) +
  tm_borders(col = "lightblue") +
  
  tm_shape(tp_vill)+tm_borders("grey")+
  
  tm_shape(tp_district) +
  tm_fill("TOWN", alpha = 0.3, title = "行政區", palette = district_colors) +   # 加透明度，不蓋住底圖
  tm_borders(lwd = 1, col = "black") +

  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("泳池500m網格熱區")

## 

## ── 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') 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_fill()`: use `fill_alpha` instead of `alpha`.
## [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>)'

免費運動場網格熱區

sport_count=as.numeric(sport_count$sport_count)
grid.sport<-grid
grid.sport$Gi_sport=localG(sport_count,grid.nb.w.in)
grid.sport$colG_sport="grey99"

#Bonferroni多重檢定修正
grid.sport$colG_sport[grid.sport$Gi_sport>=qnorm(1-0.05/n)]="purple" 

免費運動場網格熱區與行政區疊圖

tm_shape(grid.sport) +
  tm_fill("colG_sport", title = "免費運動場熱區", palette = c("grey99", "green")) +
  tm_borders(col = "lightblue") +
  
  tm_shape(tp_vill)+tm_borders("grey")+
  
  tm_shape(tp_district) +
  tm_fill("TOWN", alpha = 0.3, title = "行政區", palette = district_colors) +   # 加透明度，不蓋住底圖
  tm_borders(lwd = 1, col = "black") +

  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("免費運動場500m網格熱區")

## 

## ── 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') 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_fill()`: use `fill_alpha` instead of `alpha`.
## [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>)'

飲料店網格熱區

drink_count=as.numeric(drink_count$drink_count)
grid.drink<-grid
grid.drink$Gi_drink=localG(drink_count,grid.nb.w.in)
grid.drink$colG_drink="grey99"

#Bonferroni多重檢定修正
grid.drink$colG_drink[grid.drink$Gi_drink>=qnorm(1-0.05/n)]="orange" 

飲料店網格熱區與行政區疊圖

tm_shape(grid.drink) +
  tm_fill("colG_drink", title = "飲料店熱區", palette = c("grey99", "green")) +
  tm_borders(col = "lightblue") +
  
  tm_shape(tp_vill)+tm_borders("grey")+
  
  tm_shape(tp_district) +
  tm_fill("TOWN", alpha = 0.3, title = "行政區", palette = district_colors) +   
  tm_borders(lwd = 1, col = "black") +

  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("飲料店500m網格熱區")

## 

## ── 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') 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_fill()`: use `fill_alpha` instead of `alpha`.
## [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>)'

速食店，飲料店，免費運動場，公園與游泳池冷熱區

LISA

vill.nb=poly2nb(tp_vill)
vill.nb.w=nb2listw(vill.nb,
                 zero.policy = T)

# join fastfood_count
fastfood_joined=st_join(tp_food,tp_vill,join=st_within)

fastfood_count=fastfood_joined %>% st_drop_geometry() %>% group_by(VILLAGE) %>%
  summarise(fastfood_count=n(),.groups = "drop")

tp_vill=tp_vill %>% left_join(fastfood_count,by="VILLAGE")

tp_vill$fastfood_count[is.na(tp_vill$fastfood_count)]=0

#join drink_count

drink_joined=st_join(drink,tp_vill,join=st_within)

drink_count=drink_joined %>% st_drop_geometry() %>% group_by(VILLAGE.x) %>% 
  summarise(drink_count=n())

tp_vill=tp_vill %>% left_join(drink_count,by=c("VILLAGE"="VILLAGE.x"))

tp_vill$drink_count[is.na(tp_vill$drink_count)]=0

# join park_count
park_joined=st_join(park,tp_vill,join=st_within)

park_count=park_joined %>% st_drop_geometry() %>% group_by(VILLAGE) %>%
  summarise(park_count=n(),.groups = "drop")

tp_vill=tp_vill %>% left_join(park_count,by="VILLAGE")

tp_vill$park_count[is.na(tp_vill$park_count)]=0

# join pool_count

pool_joined=st_join(tp_pool,tp_vill,join=st_within)

pool_count=pool_joined %>% st_drop_geometry() %>% group_by(VILLAGE) %>%
  summarise(pool_count=n(),.groups = "drop")

tp_vill=tp_vill %>% left_join(pool_count,by="VILLAGE")

tp_vill$pool_count[is.na(tp_vill$pool_count)]=0

# join sport_count

sport_joined=st_join(sport,tp_vill,join=st_within)

sport_count=sport_joined %>% st_drop_geometry() %>% group_by(VILLAGE) %>%
  summarise(sport_count=n(),.groups = "drop")

tp_vill=tp_vill %>% left_join(sport_count,by="VILLAGE")

tp_vill$sport_count[is.na(tp_vill$sport_count)]=0

#village area
tp_vill$area=st_area(tp_vill)
library(units)

## udunits database from C:/Users/user/AppData/Local/R/win-library/4.4/units/share/udunits/udunits2.xml

tp_vill$area=set_units(tp_vill$area,km^2)

# density
tp_vill$fast_dens=tp_vill$fastfood_count/tp_vill$area

tp_vill$drink_dens=tp_vill$drink_count/tp_vill$area

tp_vill$park_dens=tp_vill$park_count/tp_vill$area

tp_vill$pool_dens=tp_vill$pool_count/tp_vill$area

tp_vill$sport_dens=tp_vill$sport_count/tp_vill$area

LISA_desity:每個「里」單位面積內有幾家(1/km^2)

fastfood=as.numeric(tp_vill$fast_dens)

LISA=localmoran(fastfood, vill.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=fastfood-mean(fastfood) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_vill$colI=col



tm_shape(tp_vill)+tm_polygons('colI')+tm_add_legend(type="symbol",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_borders("grey")+tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每個里的速食店密度")

## [v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.
## [v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.

drink=as.numeric(tp_vill$drink_dens)

LISA=localmoran(drink, vill.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=drink-mean(drink) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_vill$colI=col



tm_shape(tp_vill)+tm_polygons('colI')+tm_add_legend(type="symbol",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_borders("grey")+tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每個里的飲料店密度")

## [v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.
## [v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.

park=as.numeric(tp_vill$park_dens)

LISA=localmoran(park, vill.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=park-mean(park) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_vill$colI=col



tm_shape(tp_vill)+tm_polygons('colI')+tm_add_legend(type="symbol",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_borders("grey")+tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每個里的公園密度")

## [v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.
## [v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.

pool=as.numeric(tp_vill$pool_dens)

LISA=localmoran(pool, vill.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=pool-mean(pool) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_vill$colI=col



tm_shape(tp_vill)+tm_polygons('colI')+tm_add_legend(type="symbol",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_borders("grey")+tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每個里的游泳池密度")

## [v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.
## [v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.

sport=as.numeric(tp_vill$sport_dens)

LISA=localmoran(sport, vill.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=sport-mean(sport) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_vill$colI=col



tm_shape(tp_vill)+tm_polygons('colI')+tm_add_legend(type="symbol",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_borders("grey")+tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每個里免費運動場密度")

## [v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.
## [v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.

“區”為單位的冷熱區

LISA

fastfood_dens

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(fastfood_count=sum(fastfood_count))
tp_district.nb=poly2nb(tp_district)
tp_district.nb.w=nb2listw(tp_district.nb,
                          zero.policy = T)

tp_district$area=st_area(tp_district)
tp_district$area=set_units(tp_district$area,km^2)

tp_district$fast_dens=tp_district$fastfood_count/tp_district$area

fast_dens=as.numeric(tp_district$fast_dens)

LISA=localmoran(fast_dens, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=fast_dens-mean(fast_dens) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_district$colI=col

tm_shape(tp_district)+tm_polygons('colI')+tm_add_legend("fill",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_title("速食店密度熱區：大安區")

飲料店區密度冷熱區

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(drink_count=sum(drink_count))
tp_district.nb=poly2nb(tp_district)
tp_district.nb.w=nb2listw(tp_district.nb,
                          zero.policy = T)

#
tp_district$area=st_area(tp_district)
tp_district$area=set_units(tp_district$area,km^2)

tp_district$drink_dens=tp_district$drink_count/tp_district$area

drink_dens=as.numeric(tp_district$drink_dens)

LISA=localmoran(drink_dens, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=drink_dens-mean(drink_dens) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_district$colI=col

tm_shape(tp_district)+tm_polygons('colI')+tm_add_legend("fill",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_title("飲料店密度熱區_沒冷熱區")

park_dens

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(park_count=sum(park_count))
tp_district.nb=poly2nb(tp_district)
tp_district.nb.w=nb2listw(tp_district.nb,
                          zero.policy = T)

park=tp_district$park_count
LISA=localmoran(park, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

tp_district$area=st_area(tp_district)
tp_district$area=set_units(tp_district$area,km^2)

tp_district$park_dens=tp_district$park_count/tp_district$area

park_dens=as.numeric(tp_district$park_dens)

LISA=localmoran(park_dens, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=park_dens-mean(park_dens) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_district$colI=col

tm_shape(tp_district)+tm_polygons('colI')+tm_add_legend("fill",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_title("公園區密度_沒冷熱區")

pool_dens

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(pool_count=sum(pool_count))
tp_district.nb=poly2nb(tp_district)
tp_district.nb.w=nb2listw(tp_district.nb,
                          zero.policy = T)

pool=tp_district$pool_count
LISA=localmoran(pool, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

tp_district$area=st_area(tp_district)
tp_district$area=set_units(tp_district$area,km^2)

tp_district$pool_dens=tp_district$pool_count/tp_district$area

pool_dens=as.numeric(tp_district$pool_dens)

LISA=localmoran(pool_dens, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=pool_dens-mean(pool_dens) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_district$colI=col

tm_shape(tp_district)+tm_polygons('colI')+tm_add_legend("fill",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_title("游泳池區密度熱區：中正區")

sport_dens

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(sport_count=sum(sport_count))
tp_district.nb=poly2nb(tp_district)
tp_district.nb.w=nb2listw(tp_district.nb,
                          zero.policy = T)

sport=tp_district$sport_count
LISA=localmoran(sport, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

tp_district$area=st_area(tp_district)
tp_district$area=set_units(tp_district$area,km^2)

tp_district$sport_dens=tp_district$sport_count/tp_district$area

sport_dens=as.numeric(tp_district$sport_dens)

LISA=localmoran(sport_dens, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=sport_dens-mean(sport_dens) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


tp_district$colI=col

tm_shape(tp_district)+tm_polygons('colI')+tm_add_legend("fill",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_title("免費運動場區密度熱區：中正區 萬華區")

## [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.

台北市學生資料

tp_student=read.csv("C:/Users/user/Downloads/tp_students.csv")

tp_sch_stu<-as.data.frame(tp_school)
tp_sch_stu=tp_sch_stu %>% left_join(tp_student,by=c("schoolname"="學校"))

tp_sch_stu=tp_sch_stu[!is.na(tp_sch_stu$student_count),]

tp_sch_stu_sf<-st_as_sf(tp_sch_stu)

學生數分布有顯著空間自相關。

tp_district=tp_vill %>% group_by(TOWN) %>% summarise(area=sum(area))
district_stu <- st_join(tp_district, tp_sch_stu_sf, join = st_contains) %>%
  group_by(TOWN) %>%
  summarise(stu_sum = sum(student_count, na.rm = TRUE))

district_stu=left_join(district_stu, st_drop_geometry(tp_district),by="TOWN")

district_stu$stu_dens=district_stu$stu_sum/district_stu$area

district_stu=st_as_sf(district_stu)

臺北市學生數分布

vill_stu <- st_join(tp_vill, tp_sch_stu_sf, join = st_contains) %>%
  group_by(VILLAGE) %>%
  summarise(stu_sum = sum(student_count, na.rm = TRUE))

tm_shape(vill_stu) +
  tm_fill("stu_sum", palette = "OrRd", title = "每個里的學生數",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_borders("grey")+
  tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每里的學生數")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named
## "brewer.or_rd"
## Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

tp_check=st_join(tp_school, tp_vill, join = st_within)

每里學生密度

#算學生密度
vill_stu$area=st_area(vill_stu)
vill_stu$area=set_units(vill_stu$area,km^2)
vill_stu$stu_dens=vill_stu$stu_sum/vill_stu$area


tm_shape(vill_stu) +
  tm_fill("stu_dens", palette = "OrRd", title = "每個里的學生密度",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_borders("grey")+
  tm_shape(tp_district)+tm_borders("black",lwd=1.5)+tm_title("每里的學生密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named
## "brewer.or_rd"
## Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

每區的學生密度(1/km^2)

tm_shape(district_stu) +
  tm_fill("stu_dens", palette = "OrRd", title = "每區的學生密度(1/km^2)",style="jenks")+
  tm_borders(col = "black") +
  tm_layout(legend.outside = TRUE, frame = FALSE)+tm_title("每區的學生密度")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_fill()`: instead of `style = "jenks"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named
## "brewer.or_rd"
## Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

paste("區學生密度最高前三：")

## [1] "區學生密度最高前三："

paste("大同區",round(district_stu$stu_dens[district_stu$TOWN=="大同區"]))

## [1] "大同區 3234"

paste("大安區",round(district_stu$stu_dens[district_stu$TOWN=="大安區"]))

## [1] "大安區 3178"

paste("中正區",round(district_stu$stu_dens[district_stu$TOWN=="中正區"]))

## [1] "中正區 2955"

學生冷熱區

stu_dens=as.numeric(district_stu$stu_dens)

LISA=localmoran(stu_dens, tp_district.nb.w,
                zero.policy = T,
                alternative = "two.sided")

z=LISA[,4]
p=LISA[,5]
diff=stu_dens-mean(stu_dens) #自己比平均是H/L
col=c()
col[diff>0 & z>0]="red" #H/H
col[diff<0 & z>0]="blue" #L/L
col[diff>0 & z<0]="pink" #H/L
col[diff<0 & z<0]="lightblue" #L/H
col[p>0.05]="grey90" #不顯著
col[is.na(col)] = "blue"


district_stu$colI=col

tm_shape(district_stu)+tm_polygons('colI')+tm_add_legend("fill",labels=c("HH","LL","HL","LH","NS"),
col=c("red","blue","pink","lightblue","grey90"))+tm_title("學生密度冷熱區_萬華區(L-H)")