開発経済学とその他応用分野を学ぶ院生

人間万事塞翁が馬を大切にしている応用経済学徒. 2020年4月から開発・計量・プログラミング関連の記事を書きます.

ggplot2::geom_pointで人口に関する地図を書く(その1)

はじめに

この記事では、以下のような関西地域の人口に関する地図をRで作成する方法を紹介する。海外新聞ぽいフォントとかを使いお洒落に仕上げることが目標。出来るだけ丁寧に書いているつもり。全てのコードを載せている。ステップバイステップでやるのを心がけた。バリエーションがないため、レンジを区切って離散クラスにすればもっと綺麗に描画できるが、疲れたので次回以降にする。(間違いなどあればコメントください。)

 

f:id:econgrad:20201004160901p:plain

 

 

1. データソース

国土数値情報のウェブサイトから「1kmメッシュ別将来推計人口(H30国政局推計)(shape形式版)」の関西二府四県のデータをダウンロードする。ファイルサイズが大きいので注意。
 

f:id:econgrad:20201004011848p:plain
 

ダウンロードしたshapeファイルをQGIS上で開くとこんな感じになる(EPSG:6674)。    

f:id:econgrad:20201004012951p:plain
 

 

2. データ読み込み & 確認

# Setup
library(ggplot2)
library(dplyr)
library(sf)
# library(ggthemes)
library(viridis)
library(extrafont)

# Set your wd
setwd()

# Load the data
shiga <- read_sf("1km_mesh_suikei_2018_shape_25/1km_mesh_2018_25.shp")
kyoto <- read_sf("1km_mesh_suikei_2018_shape_26/1km_mesh_2018_26.shp")
osaka <- read_sf("1km_mesh_suikei_2018_shape_27/1km_mesh_2018_27.shp")
hyogo <- read_sf("1km_mesh_suikei_2018_shape_28/1km_mesh_2018_28.shp")
nara <- read_sf("1km_mesh_suikei_2018_shape_29/1km_mesh_2018_29.shp")
wakayama <- read_sf("1km_mesh_suikei_2018_shape_30/1km_mesh_2018_30.shp")

 

試しに先ほどの図と同じようなものを書いてみる。

 

# Draw Maps
kansai <- ggplot() + 
  geom_sf(data = shiga, color = "red") +
  geom_sf(data = kyoto, color = "green") +
  geom_sf(data = osaka, color = "blue") +
  geom_sf(data = hyogo, color = "cyan") +
  geom_sf(data = nara, color = "magenta") +
  geom_sf(data = wakayama, color = "yellow") +
  theme_minimal() + theme(panel.grid = element_blank()) + 
  labs(x = "Longitude", y = "Latitude") +
  theme(axis.text = element_text(size = 10),
        axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10))
kansai
ggsave(kansai, filename = "kansai.png", width = 6, height = 8)

 

こんな感じになる。色は適当。後で描画するものを映えさせるためにあえてダサくした。

 

f:id:econgrad:20201004160713p:plain

 

3. データの変換etc

上で見たように、今読み込んでいるデータは1kmメッシュのデータである。これをポイント(点)に変換したい。重心を取る。
sf::st_centroidを用いる。

 

# Calculate centroid of polygon using sf::st_centroid
shiga$centroids <- st_transform(shiga, 6674) %>% 
  st_centroid() %>% 
  st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>%
  st_geometry() 

kyoto$centroids <- st_transform(kyoto, 6674) %>% 
  st_centroid() %>% 
  st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>%
  st_geometry() 

osaka$centroids <- st_transform(osaka, 6674) %>% 
  st_centroid() %>% 
  st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>%
  st_geometry() 

hyogo$centroids <- st_transform(hyogo, 6674) %>% 
  st_centroid() %>% 
  st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>%
  st_geometry() 

nara$centroids <- st_transform(nara, 6674) %>% 
  st_centroid() %>% 
  st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>%
  st_geometry() 

wakayama$centroids <- st_transform(wakayama, 6674) %>% 
  st_centroid() %>% 
  st_transform(., '+proj=longlat +ellps=GRS80 +no_defs') %>%
  st_geometry() 

 

データフレームに変換する。

 

# Convert to data frame
shiga <- as.data.frame(shiga, xy = T)
kyoto <- as.data.frame(kyoto, xy = T)
osaka <- as.data.frame(osaka, xy = T)
hyogo <- as.data.frame(hyogo, xy = T)
nara <- as.data.frame(nara, xy = T)
wakayama <- as.data.frame(wakayama, xy = T)

 

longitudeとlatitudeを抽出し、それぞれを独立の列に変数として加える。

 

# Extract and split lat long coordinates from the point
# refer: https://stackoverflow.com/questions/38637070/extract-and-split-lat-long-coordinates-from-wkt-point-data-in-r
shiga$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", shiga$centroids))
shiga$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", shiga$centroids))

kyoto$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", kyoto$centroids))
kyoto$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", kyoto$centroids))

osaka$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", osaka$centroids))
osaka$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", osaka$centroids))

hyogo$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", hyogo$centroids))
hyogo$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", hyogo$centroids))

nara$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", nara$centroids))
nara$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", nara$centroids))

wakayama$long <- as.numeric(gsub(".*?([-]*[0-9]+[.][0-9]+).*", "\\1", wakayama$centroids))
wakayama$lat <- as.numeric(gsub(".* ([-]*[0-9]+[.][0-9]+).*", "\\1", wakayama$centroids))

 

必要のない変数を削除する。everything()で順番を変更する。これは好みの問題。Stataのorderコマンドと同じ。

 

shiga <- shiga %>%
  dplyr::select(long, lat, PTN_2015) %>%
  dplyr::select(long, lat, PTN_2015, everything()) %>%
  dplyr::rename(Population = PTN_2015)

kyoto <- kyoto %>%
  dplyr::select(long, lat, PTN_2015) %>%
  dplyr::select(long, lat, PTN_2015, everything()) %>%
  dplyr::rename(Population = PTN_2015)

osaka <- osaka %>%
  dplyr::select(long, lat, PTN_2015) %>%
  dplyr::select(long, lat, PTN_2015, everything()) %>%
  dplyr::rename(Population = PTN_2015)

hyogo <- hyogo %>%
  dplyr::select(long, lat, PTN_2015) %>%
  dplyr::select(long, lat, PTN_2015, everything()) %>%
  dplyr::rename(Population = PTN_2015)

nara <- nara %>%
  dplyr::select(long, lat, PTN_2015) %>%
  dplyr::select(long, lat, PTN_2015, everything()) %>%
  dplyr::rename(Population = PTN_2015)

wakayama <- wakayama %>%
  dplyr::select(long, lat, PTN_2015) %>%
  dplyr::select(long, lat, PTN_2015, everything()) %>%
  dplyr::rename(Population = PTN_2015)

# Append the prefectures' data
kansai_pop <- dplyr::bind_rows(
  shiga, kyoto, osaka, hyogo, nara, wakayama
)

 

4. Maps−地図を書く−

冒頭で掲載した地図を作成する。

 

# Draw a map (1)
pop_map1 <- ggplot() + 
  geom_point(data= kansai_pop, 
             aes(x = long, y = lat, color = Population),
             alpha = 0.5) + 
  scale_color_gradient(low = "#edd6d4", high = "#7A0018") + 
  labs(x = NULL, 
       y = NULL, 
       title = "Kansai's Regional Demographics", 
       subtitle = "Population in Kansai Municipalities, 2015", 
       caption = "Source: National Spatial Planning and Regional Policy Bureau") +
  theme_void() +
  theme(
    plot.title = element_text(vjust = -2, hjust = 0.1),
    plot.subtitle = element_text(vjust = -2, hjust = 0.1),
    plot.caption = element_text(vjust = 3, hjust = 0.8),
    text = element_text(family = "Optima"),
    plot.background = element_rect(fill = "#f5f5f2"), 
    legend.position = c(0.13, 0.18))

pop_map1
ggsave(pop_map1, filename = "pop_map1.png", width = 6, height = 8)

 

f:id:econgrad:20201004160901p:plain

 

せっかくgeom_pointを使っているので、ポイントのサイズを人口に合わせて変えてみる。

 

# Draw a map (1)
pop_map2 <- ggplot() + 
  geom_point(data= kansai_pop, 
             aes(x = long, y = lat,
                 size = Population, color = Population),
             alpha = 0.3) + 
  scale_color_gradient(low = "#edd6d4", high = "#7A0018") + 
  # Show a legend for a color
  guides(size = FALSE) + 
  scale_size_area(max_size = 15) + 
  labs(x = NULL, 
       y = NULL, 
       title = "Kansai's Regional Demographics", 
       subtitle = "Population in Kansai Municipalities, 2015", 
       caption = "Source: National Spatial Planning and Regional Policy Bureau") +
  theme_void() +
  theme(
    plot.title = element_text(vjust = -2, hjust = 0.1),
    plot.subtitle = element_text(vjust = -2, hjust = 0.1),
    plot.caption = element_text(vjust = 3, hjust = 0.8),
    text = element_text(family = "Optima"),
    plot.background = element_rect(fill = "#f5f5f2"), 
    legend.position = c(0.13, 0.18))

pop_map2
ggsave(pop_map2, filename = "pop_map2.png", width = 6, height = 8)

 

f:id:econgrad:20201004163438p:plain

 

参考にしたサイト timogrossenbacher.ch

econgrad.hatenablog.com

econgrad.hatenablog.com

econgrad.hatenablog.com

econgrad.hatenablog.com

econgrad.hatenablog.com

econgrad.hatenablog.com