日立市の現状

Package・ディレクトリの準備

library(sf)
library(tidyverse)

setwd("~/R/ForTeaching")

データ読み込み

国土数値情報　>>　500mメッシュ別将来推計人口（H30国政局推計）（shape形式版）　>>　茨城県　>>　世界測地系

#行政区域（茨城県）
Ibaraki_map<-
  read_sf("Hitachi/N03-20210101_08_GML/N03-21_08_210101.shp")

#日立市選別
Ibaraki_map %>%
  filter(N03_004=="日立市") -> 
  Hitachi_shi_map 

#茨城メッシュ
Ibaraki_mesh<-
  read_sf("Hitachi/500m_mesh_suikei_2018_shape_08/500m_mesh_2018_08.shp") 

#日立市選別
Ibaraki_mesh %>%
  filter(SHICODE==8202) -> 
  Hitachi_shi_mesh

人口データの可視化

500mメッシュデータで可視化する。

#500mメッシュ別人口の可視化
ggplot()+
  geom_sf(data=Hitachi_shi_map, fill="white")+ 
  geom_sf(data=Hitachi_shi_mesh, 
          aes(fill=PTN_2015))+
  scale_fill_distiller(palette="YlGnBu", direction=1)+
  geom_sf(data=Hitachi_shi_map, fill="NA")+ 
  labs(fill="人口", 
       caption="出典：国土交通省国土数値情報")+
  ggtitle("日立市500mメッシュ別人口（2015）")+
  theme_gray (base_family = "HiraKakuPro-W3")

# 日本語保存のやり方
# 日本語フォントの扱いで
# theme_gray (base_family = "HiraKakuPro-W3")
# を設定するとPDF保存できない。しかし、
# theme_gray (base_family = "HiraKakuPro-W3")
# を設定しないとhtml表示で文字化けする。
# そのため以下では
# theme_gray (base_family = "HiraKakuPro-W3")
# だけ外して、PDFとしてggsaveで保存するためだけに再度生成している

ggplot()+
  geom_sf(data=Hitachi_shi_map, fill="white")+ 
  geom_sf(data=Hitachi_shi_mesh, 
          aes(fill=PTN_2015))+
  scale_fill_distiller(palette="YlGnBu", direction=1)+
  geom_sf(data=Hitachi_shi_map, fill="NA")+ 
  labs(fill="人口", 
       caption="出典：国土交通省国土数値情報")+
  ggtitle("日立市500mメッシュ別人口（2015）")->HitachiPop500meshPDF

ggsave(filename="Hitachi/HitachiPop500mesh.pdf",
       device=cairo_pdf,
       family="Hiragino Sans",
       plot=HitachiPop500meshPDF,
       unit="in",width=6,height=10)

森林エリアの描写

国土数値情報 >> 1. 国土の「土地利用」 >> 森林地域（ポリゴン）>> 茨城県 (H27)

をlocalに保存し、ForestAreaIbarakiというフォルダ名の変更した。

ForestIbaraki<-read_sf("ForestAreaIbaraki/a001080020160207.shp")

ggplot()+
  geom_sf(data=ForestIbaraki,fill="green")

日立市を取り出す

# ポリゴンの頂点が重複していると新しいsfではエラーが出るので、
# s2をFALSEにしている
# https://stackoverflow.com/questions/68808238/how-to-fix-spherical-geometry-errors-caused-by-conversion-from-geos-to-s2
sf_use_s2(FALSE) 


# JGD2000で投影法を定義 (JGD2000に対応するEPSGは4612)

Hitachi_shi_map<-Hitachi_shi_map %>%
  sf::st_set_crs(sf::st_crs(4612))

ForestIbaraki<-ForestIbaraki %>%
  sf::st_set_crs(sf::st_crs(4612))

ForestHitachi<-st_sf(st_intersection(st_union(Hitachi_shi_map),
                                     st_union(ForestIbaraki)))

ggplot()+
  geom_sf(data=ForestHitachi,fill="green")

最後に森林と人口を合わせる。人が住んでいないエリアは森林であることがわかる。

ggplot()+
  geom_sf(data=Hitachi_shi_map, fill="white")+ 
  geom_sf(data=ForestHitachi, fill="forestgreen")+
  geom_sf(data=Hitachi_shi_mesh, 
          aes(fill=PTN_2015))+
  scale_fill_distiller(palette="YlGnBu", direction=1)+
  #geom_sf(data=Hitachi_shi_map, fill="NA")+ 
  labs(fill="人口", 
       caption="出典：国土交通省国土数値情報")+
  ggtitle("日立市500mメッシュ別人口（2015）")+
  theme_gray (base_family = "HiraKakuPro-W3")

急傾斜地崩壊危険区域の追加

国土数値情報 >> 2. 政策区域 >> 災害・防災 >> 茨城県 (R2)

をlocalに保存し、Slopeというフォルダ名の変更した。

急傾斜地崩壊危険区域の指定を要する土地（区域）は、以下の[1]及び[2]の区域を包括する区域である(出所：国土交通省 )。

• [1]崩壊するおそれのある急傾斜地（傾斜度が30度以上の土地をいう。以下同じ。）で、その崩壊により相当数の居住者その他の者に被害のおそれのあるもの

• [2]に隣接する土地のうち、急傾斜地の崩壊が助長・誘発されるおそれがないようにするため、一定の行為制限の必要がある土地の区域

DangerSlopeIbaraki<-read_sf("Slope/A47-21_08.shp",
                            options = c("encoding=CP932"))

ggplot()+
  geom_sf(data=Ibaraki_map)+
  geom_sf(data=DangerSlopeIbaraki,color="red")

DangerSlopeHitachi<-DangerSlopeIbaraki |>
  filter(A47_003=="日立市")

日立市の地図に追加。これをみると急傾斜地崩壊危険区域が人口が多いエリアにもあることから、森林エリア以外にも傾斜のある地形が多いことがわかる。

ggplot()+
  geom_sf(data=Hitachi_shi_map, fill="white")+ 
  geom_sf(data=ForestHitachi, fill="forestgreen")+
  geom_sf(data=Hitachi_shi_mesh, 
          aes(fill=PTN_2015))+
  scale_fill_distiller(palette="YlGnBu", direction=1)+
  geom_sf(data=DangerSlopeHitachi, color="red")+ 
  labs(fill="人口", 
       caption="出典：国土交通省国土数値情報\n (赤点は急傾斜地崩壊危険区域)")+
  ggtitle("日立市500mメッシュ別人口（2015）")+
  theme_gray (base_family = "HiraKakuPro-W3")

急傾斜地崩壊危険区域の人口を集計する

st_intersection()は、形を切り抜くが、st_intersects()は、重なっているかどうかを論理で返すリストを生成する。

# 投影法をあわせる
Hitachi_shi_mesh<-Hitachi_shi_mesh %>%
  sf::st_set_crs(sf::st_crs(4612))

# 投影法をあわせる
DangerSlopeHitachi<-DangerSlopeHitachi %>%
  sf::st_set_crs(sf::st_crs(4612))

# 2つのポリゴンの交わるところを抽出
#DangerPopHitachi<-st_intersection(Hitachi_shi_mesh,
#                                  DangerSlopeHitachi)
#

# 2つのポリゴンの交わるところを抽出
DangerPopHitachi<-st_intersects(Hitachi_shi_mesh,
                                  DangerSlopeHitachi)

# 結果のリストで重なりがあったIDを抽出
DPH_logical = lengths(DangerPopHitachi) > 0
# メッシュのうち、重なりがあった部分を抽出
Hitachi_shi_mesh2<-Hitachi_shi_mesh[DPH_logical,]

# 上記と同じ結果を得る別の方法
# https://geocompr.robinlovelace.net/spatial-operations.html
Hitachi_shi_mesh3 <-Hitachi_shi_mesh |>
  st_filter(y=DangerSlopeHitachi, .predicate=st_intersects)

ggplot()+
  geom_sf(data=Hitachi_shi_map, fill="white")+ 
  #geom_sf(data=ForestHitachi, fill="forestgreen")+
  geom_sf(data=Hitachi_shi_mesh2, 
          aes(fill=PTN_2015))+
  scale_fill_distiller(palette="YlGnBu", direction=1)+
  #geom_sf(data=DangerSlopeHitachi, color="red")+ 
  labs(fill="人口", 
       caption="出典：国土交通省国土数値情報\n (赤点は急傾斜地崩壊危険区域)")+
  ggtitle("急傾斜地崩壊危険区域と重なる\n日立市500mメッシュ別人口（2015）")+
  theme_gray (base_family = "HiraKakuPro-W3")

日立市では市の人口の約10分の1が急傾斜地崩壊危険区域に住んでいる。

# 人数
DangerPopHitachiSum15<-sum(Hitachi_shi_mesh2$PTN_2015)
DangerPopHitachiSum20<-sum(Hitachi_shi_mesh2$PTN_2020)
DangerPopHitachiSum25<-sum(Hitachi_shi_mesh2$PTN_2025)
DangerPopHitachiSum30<-sum(Hitachi_shi_mesh2$PTN_2030)
DangerPopHitachiSum35<-sum(Hitachi_shi_mesh2$PTN_2035)
DangerPopHitachiSum40<-sum(Hitachi_shi_mesh2$PTN_2040)

Year<-c(2015,2020,2025,2030,2035,2040)
DP<-c(DangerPopHitachiSum15,DangerPopHitachiSum20,
      DangerPopHitachiSum25,DangerPopHitachiSum30,
      DangerPopHitachiSum35,DangerPopHitachiSum40)

DPHitachiSum<-as.data.frame(cbind(Year, DP))

ggplot()+
  geom_col(data=DPHitachiSum,aes(x=Year,y=DP),fill="firebrick3")+
  xlab("年")+ylab("急傾斜地崩壊危険区域内の人口 (単位:人)")+
    labs(caption="出典：国土交通省国土数値情報")+
  theme_gray (base_family = "HiraKakuPro-W3")

# 市人口に占める割合
DPHitachiRatio15<-DangerPopHitachiSum15/sum(Hitachi_shi_mesh$PTN_2015)
DPHitachiRatio20<-DangerPopHitachiSum20/sum(Hitachi_shi_mesh$PTN_2020)
DPHitachiRatio25<-DangerPopHitachiSum25/sum(Hitachi_shi_mesh$PTN_2025)
DPHitachiRatio30<-DangerPopHitachiSum30/sum(Hitachi_shi_mesh$PTN_2030)
DPHitachiRatio35<-DangerPopHitachiSum35/sum(Hitachi_shi_mesh$PTN_2035)
DPHitachiRatio40<-DangerPopHitachiSum40/sum(Hitachi_shi_mesh$PTN_2040)

DPR<-c(DPHitachiRatio15,DPHitachiRatio20,
      DPHitachiRatio25,DPHitachiRatio30,
      DPHitachiRatio35,DPHitachiRatio40)

DPHitachiRatio<-as.data.frame(cbind(Year, DPR))

ggplot()+
  geom_line(data=DPHitachiRatio,aes(x=Year,y=DPR),
            color="firebrick3",size=1)+
  geom_point(data=DPHitachiRatio,aes(x=Year,y=DPR),
             color="firebrick3",size=2)+
  ylim(0,0.15)+
  xlab("年")+ylab("急傾斜地崩壊危険区域内の人口割合")+
  labs(caption="出典：国土交通省国土数値情報")+
  theme_gray (base_family = "HiraKakuPro-W3")