目的

土地利用細分メッシュデータを用いて土地利用状況を可視化.土地利用の現況が衛星画像を用いて土地分類基準に従い判読されている.

  • 100mメッシュのため,観測数が多く(例えば佐賀市の地図の作成のために117万メッシュを取り扱う),地図にする際には時間がかかることに注意.
  • ここでは取り組まない(時間かかるため)が,時系列に地図を並べることより,土地利用がどのように変化したかを知ることができる.

ライブラリ:sftidyverse

library(sf)
library(tidyverse)

データ整理

シェープファイルのダウンロード

最初に,佐賀市の行政区域をダウンロード.

  • 国土数値情報ダウンロードサービス > 行政区域(ポリゴン) > 佐賀県 > 世界測地系 令和4年

    • ダウンロードしたZIPファイルは同じフォルダ(ディレクトリ)に納め,ZIPファイルを選択 > 右クリック > すべて展開 > 展開先の選択とファイルの選択 > 展開を実行する.展開後は最初にダウンロードしたZIPファイルを削除.
    • シェープファイルの読込(read_sf()):行政区域をSaga_mapと名付ける.
      • Saga_mapから佐賀市を抜き出す(filter()).
#行政区域(佐賀県)
Saga_map <-
  read_sf("N03-20220101_41_GML/N03-22_41_220101.shp")

#佐賀市選別
Saga_map %>%
  filter(N03_004=="佐賀市") -> 
  Saga_shi_map 

続いて,佐賀市周辺の土地利用細分メッシュデータをダウンロード.

Saga_tochiriyo_map1 <- 
  read_sf("L03-b-21_4930-jgd2011_GML/L03-b-21_4930.shp")

Saga_tochiriyo_map2 <- 
  read_sf("L03-b-21_5030-jgd2011_GML/L03-b-21_5030.shp")

データの中身の確認

head()を用いてデータの中身を確認.

  • L03b_001はメッシュコードを,L03b_002土地利用種別コードを,L03b_003は衛星写真撮影年月日を示す.
Saga_tochiriyo_map1 %>% 
  head()
L03b_001 L03b_002 L03b_003 geometry
4930010000 1500 20201104 POLYGON ((130.125 32.66667,…
4930010001 1500 20201104 POLYGON ((130.1263 32.66667…
4930010002 1500 20201104 POLYGON ((130.1275 32.66667…
4930010003 0500 20201104 POLYGON ((130.1288 32.66667…
4930010004 0200 20201104 POLYGON ((130.13 32.66667, …
4930010005 0200 20201104 POLYGON ((130.1313 32.66667…

土地利用種別

土地利用種別を地図上に色分けする.ここでは,次のように8分類に指定(色は好み).

土地利用種別
土地利用種別 コード
田・その他農用地 0100,0200 #DEF5E5FF
森林 0500 #ADE3C0FF
荒れ地 0600 #6CD3ADFF
建物用地 0700 #366DA0FF
幹線交通用地 0901, 0902 #7C7B78FF
その他の用地 1000 #3E356BFF
河川池・湖沼・海浜 1100,1400,1500 #35264CFF
ゴルフ場 1600 #AADC32FF

凡例準備

シェープファイルには土地利用種別がコードで分類されている.このままだと,凡例もコードで分類されるため,そのコードが意味することをすぐには判断できない.そこで,コードの意味を示す列を作成.

  • 最初に,土地利用種別とコードの対応関係を指示(my_listと名付ける).
  • 次に新しい列「土地利用種別」を作成する.fct_recode()を用いて,L03b_002列にある元のコードを新しい土地利用種別に変更.
    • ここでは,my_listをすでに作成しているため,これを利用(!!!my_list).
    • よく使うR小技集(アクセス日:2022-08-12;再アクセス日:2024-03-08)参照.
#土地利用種別とコードの対応関係
my_list<-
  c("田・その他農用地"="0100",
    "田・その他農用地"="0200",
    "森林"="0500",
    "荒れ地"="0600",
    "建物用地"="0700",
    "道路・鉄道"="0901",
    "道路・鉄道"="0902",
    "その他の用地"="1000",
    "河川池・湖沼・海浜"="1100",
    "河川池・湖沼・海浜"="1400",
    "河川池・湖沼・海浜"="1500",
    "ゴルフ場"="1600")

Saga_tochiriyo_map1 %>%
  mutate(土地利用種別=fct_recode(L03b_002, !!!my_list))->
  Saga_tochiriyo_map1

Saga_tochiriyo_map2 %>% 
  mutate(土地利用種別=fct_recode(L03b_002, !!!my_list))->
  Saga_tochiriyo_map2

head()を用いて土地利用種別を示す列ができているか確認.

Saga_tochiriyo_map1 %>%
  head()
L03b_001 L03b_002 L03b_003 geometry 土地利用種別
4930010000 1500 20201104 POLYGON ((130.125 32.66667,… 河川池・湖沼・海浜
4930010001 1500 20201104 POLYGON ((130.1263 32.66667… 河川池・湖沼・海浜
4930010002 1500 20201104 POLYGON ((130.1275 32.66667… 河川池・湖沼・海浜
4930010003 0500 20201104 POLYGON ((130.1288 32.66667… 森林
4930010004 0200 20201104 POLYGON ((130.13 32.66667, … 田・その他農用地
4930010005 0200 20201104 POLYGON ((130.1313 32.66667… 田・その他農用地

色の指定

色については上の表に合わせて,次のようにmy_color作成.

my_color=c("#DEF5E5FF",
           "#ADE3C0FF",
           "#6CD3ADFF",
           "#366DA0FF",
           "#7C7B78FF",
           "#3E356BFF",
           "#35264CFF",
           "#AADC32FF")

完成図(佐賀市)

3枚の地図を重ねる.

  • 1,2枚目.土地利用種別.
    • 土地利用種別にメッシュの色分けする(aes()).メッシュの枠,枠内を同色するため,color(枠)とfill(枠内)を指示.
    • scale_color_manual()で枠の色をマニュアルで指定.色はmy_colorを利用.guide="legend"で凡例を整える.
    • scale_fill_manual()で枠内の色をマニュアルで指定.色はmy_colorを利用.
  • 3枚目.佐賀市の枠を可視化.
    • 枠の色は#FCFDBFFFを指示.線の太さ(linewidth)は0.8.枠内は無色(NA
  • theme(legend.position="bottom")で凡例を地図の下に配置.
ggplot()+
  geom_sf(data=Saga_tochiriyo_map1,
          aes(color=土地利用種別,
              fill=土地利用種別),
          size=0.01)+
  geom_sf(data=Saga_tochiriyo_map2,
          aes(color=土地利用種別,
              fill=土地利用種別),
          size=0.01)+
  scale_color_manual(values=my_color,
                     guide="legend")+ 
  scale_fill_manual(values=my_color)+
  geom_sf(data=Saga_shi_map, fill="NA", 
          color="#FCFDBFFF", linewidth=0.8)+
  labs(caption="出典:国土交通省国土数値情報")+
  ggtitle("佐賀市周辺の土地利用状況(2021年)")+
  coord_sf(xlim= c(130.1, 130.6), ylim=c(33.1, 33.4))+
  theme_bw()+
  theme(legend.position="bottom")

弘前市

同様にして,弘前市周辺の土地利用状況を可視化.

  • 弘前市は最新(2022-08-12現在)が2016年(平成28年).
#行政区域(青森県,2016年)
Aomori_map <-
  read_sf("N03-160101_02_GML/N03-16_02_160101.shp")

#弘前市選別
Aomori_map %>%
  filter(N03_004=="弘前市") -> 
  Hirosaki_shi_map
#列が日本語
#オプションをつけないと文字化け
Hirosaki_tochiriyo_map1 <- 
  read_sf("L03-b-16_6040-jgd_GML/L03-b-16_6040.shp",
          options=c("encoding=CP932"))

Hirosaki_tochiriyo_map2 <- 
  read_sf("L03-b-16_6140-jgd_GML/L03-b-16_6140.shp",
          options=c("encoding=CP932"))
#土地利用種別の列作成
Hirosaki_tochiriyo_map1 %>% 
 mutate(土地利用種別=fct_recode(土地利用種, !!!my_list))->
 Hirosaki_tochiriyo_map1

Hirosaki_tochiriyo_map2 %>% 
  mutate(土地利用種別=fct_recode(土地利用種, !!!my_list))->
  Hirosaki_tochiriyo_map2

完成図(弘前市)

ggplot()+
  geom_sf(data=Hirosaki_tochiriyo_map1,
          aes(color=土地利用種別,
              fill=土地利用種別),
          size=0.01)+
  geom_sf(data=Hirosaki_tochiriyo_map2,
          aes(color=土地利用種別,
              fill=土地利用種別),
          size=0.01)+
  scale_color_manual(values=my_color,
                     guide="legend")+ 
  scale_fill_manual(values=my_color)+
  geom_sf(data=Hirosaki_shi_map, fill="NA", 
          color="#FCFDBFFF", linewidth=0.8)+
  labs(caption="出典:国土交通省国土数値情報")+
  ggtitle("弘前市周辺の土地利用状況(2016年)")+
  coord_sf(xlim= c(140.2, 140.7), ylim=c(40.5, 40.8))+
  theme_bw()+
  theme(legend.position="bottom")

Rによる地理空間データの可視化