library(tidyverse)
library(tidygeocoder)
library(sf)
library(jpgrid)
9: 地域メッシュ統計
はじめに
今回も、これまでの講義内容の復習を兼ねて、R
を使った立地分析の演習を行います。 今日の演習では、いろいろな空間データをつかって、佐賀市内のコンビニエンスストアの立地について考えてみることにします。
今日の分析では、新たにjpgrid
パッケージを使いますので、インストールしてください。 そして、以下のlibrary
関数を実行して、使用するパッケージをロードしてください。
コンビニエンスストアのデータ
コンビニエンスストア(以下、コンビニと呼びます)のポイントデータは、無料で手に入るものはありませんので、購入するか自分で作成するかする必要があります。 今回の演習では、佐賀市内のコンビニの住所データを用意1しましたので、これをジオコードして使いましょう。
CSVファイルをダウンロードし、プロジェクトのdataフォルダに移動してください。 これをreadr::read_csv
関数で読み込みます。
<- read_csv('data/saga_cvs.csv') dat
Rows: 110 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): brand, name, address
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(dat)
Rows: 110
Columns: 3
$ brand <chr> "ファミリーマート", "ファミリーマート", "ファミリーマート", "…
$ name <chr> "ファミリーマート佐賀駅北店", "ファミリーマート佐賀県庁店", "…
$ address <chr> "佐賀県佐賀市駅前中央2-2593", "佐賀県佐賀市城内一丁目1番59号地…
このデータには、佐賀市内の110件のコンビニの店名と位置情報が入っています。
チェーンごとのコンビニ立地件数を確認してみましょう。
|>
dat count(brand) |>
arrange(desc(n))
# A tibble: 6 × 2
brand n
<chr> <int>
1 セブンイレブン 60
2 ローソン 27
3 ファミリーマート 17
4 ミニストップ 4
5 デイリ-ヤマザキ 1
6 ポプラ 1
佐賀市内のコンビニは、セブンイレブンが最も多く、次いでローソン、ファミリーマートの順になっています。
tidygeocoder::geocode
関数を使って、ジオコード(緯度経度情報を付与)します(110件のデータがありますので、少し時間がかかります)。 そして、sf::st_as_sf
関数で、Simple Feature(ポイントデータ)に変換しましょう。 CRSはWGS84として読み込んだのち、sf::st_transform
関数で平面直角座標系II系に投影変換しましょう。
<- dat |>
cvs geocode(address, method = 'arcgis') |>
st_as_sf(coords = c('long', 'lat'), crs = 4326) |>
st_transform(6670)
作成したデータを地図にしてみましょう。
佐賀市のコンビニは、そのほとんどが市南部の平野部に集中していることがわかります2。
何らかの理由でジオコードがうまくいかなかった場合は、下のリンクからGeoJSONファイルをダウンロードして、sf::st_read
関数で読み込んでください。
<- st_read('data/saga_cvs.geojson') |>
cvs st_transform(6670)
背景地図として、佐賀市のポリゴンデータを用意します。
<-
saga_city st_read('data/410004saga.geojson') |>
filter(GST_NAME == '佐賀市') |>
st_transform(6670)
Reading layer `410004saga' from data source
`/Users/kzktmr/Documents/LectureNotes/data/410004saga.geojson'
using driver `GeoJSON'
Simple feature collection with 20 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 129.7368 ymin: 32.95054 xmax: 130.5424 ymax: 33.6189
Geodetic CRS: WGS 84
ggplot
で地図にしてみましょう。
ggplot() +
geom_sf(data = saga_city) +
geom_sf(aes(color = brand), data = cvs) +
theme_minimal()
コンビニエンスストアの周辺人口の確認
都市部のコンビニの商圏は、徒歩5分の範囲であるといわれます。 佐賀市のコンビニストアについて、商圏人口がどのくらいいるのか調べてみましょう。
人口データは国勢調査の基本単位区データを使いましょう。e-Statページにアクセスして、 基本単位区の境界データをダウンロードしてください。
- 境界データダウンロード
- [小地域]→[国勢調査]→[2020年]→[小地域(基本単位区)(JGD2011)]→[世界測地系平面直角座標系・Shapefile]→[41 佐賀県]→[41201 佐賀市]からダウンロード
- 境界データ(B002005212020XYSWC41201-JGD2011.zip)を展開してできたフォルダを、プロジェクトのdataフォルダに移動
- [小地域]→[国勢調査]→[2020年]→[小地域(基本単位区)(JGD2011)]→[世界測地系平面直角座標系・Shapefile]→[41 佐賀県]→[41201 佐賀市]からダウンロード
境界データをsf::st_read
関数を使って読み込みます。
<-
district st_read('data/B002005212020XYSWC41201-JGD2011/r2kb41201.shp')
Reading layer `r2kb41201' from data source
`/Users/kzktmr/Documents/LectureNotes/data/B002005212020XYSWC41201-JGD2011/r2kb41201.shp'
using driver `ESRI Shapefile'
Simple feature collection with 3835 features and 37 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -80084.88 ymin: 15840.25 xmax: -57844.26 ymax: 53682.39
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
コンビニの商圏人口を調べるために、
- コンビニの徒歩5分(≒350 m)のバッファーを作成
- バッファーに代表点が含まれる基本単位区の人口を集計する
という操作を行います。 まず、基本単位区はポリゴンデータですから、sf::st_centroid
関数を使って代表点のポイントデータに変換します。
<- district |> st_centroid() pop
Warning: st_centroid assumes attributes are constant over geometries
コンビニの350 mバッファーを作成します。
<- cvs |>
buffer st_buffer(dist = 350)
コンビニのバッファーと基本単位区の代表点を重ねて地図にしてみましょう。
ggplot() +
geom_sf(data = saga_city) +
geom_sf(aes(fill = JINKO), data = district, color = NA) +
geom_sf(data = pop, size = 0.1, alpha = 0.5) +
geom_sf(data = buffer, color = 'red', fill = NA) +
scale_fill_distiller(direction = 1) +
theme_minimal()
佐賀市中心部を拡大してみます。
ggplot() +
geom_sf(data = saga_city) +
geom_sf(aes(fill = JINKO), data = district, color = NA) +
geom_sf(data = pop, size = 0.1, alpha = 0.5) +
geom_sf(data = buffer, color = 'red', fill = NA) +
scale_fill_distiller(direction = 1) +
xlim(-70000, -60000) + ylim(25000, 35000) +
theme_minimal()
それぞれの店舗の350 mバッファー(商圏)に、たくさんの基本単位区(の図形重心)が含まれていることがわかると思います。
xlim
やylim
関数には、それぞれx軸とy軸の範囲(下限値と上限値)を与える必要があります。投影変換後の座標のあたりをつけるためには、st_bbox
関数が便利かもしれません。
st_bbox(saga_city)
xmin ymin xmax ymax
-80084.89 15840.25 -57844.26 53682.39
この結果と図を見比べると、およその見当がつけられると思います。
それでは、店舗の商圏に含まれる基本単位区の人口(商圏人口)を集計してみましょう。 sf::st_join
関数と、group_by
とsummarise
関数を組み合わせて使います。
sf::st_join
は、空間結合のための関数で、空間データどうしのleft_join
関数に相当します。 ここでは、コンビニ商圏のポリゴンデータに対して、そのポリゴンと交わる基本単位区代表点のポイントデータを全て結合します。 そしてgroup_by
関数を使ってコンビニごとにグルーピングした上で、summarise
関数を使って、商圏ごとに、交わる基本単位区の人口を合計しています。
<- buffer |>
market st_join(pop) |>
group_by(name) |> summarise(population = sum(JINKO))
|> arrange(desc(population)) market
Simple feature collection with 110 features and 2 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -71976.77 ymin: 21643.1 xmax: -59447.42 ymax: 47698.05
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
# A tibble: 110 × 3
name population geometry
<chr> <int> <POLYGON [m]>
1 ファミリーマートJR佐賀駅店 2600 ((-65100.58 29509.74, -65101.0…
2 ローソン佐賀若楠三丁目店 2557 ((-66230.2 30930.74, -66230.68…
3 セブンイレブン佐賀兵庫北店 2490 ((-63758.9 29986.41, -63759.38…
4 ファミリーマート佐賀駅南口店 2433 ((-65070.87 29388.87, -65071.3…
5 ローソン佐賀高伝寺前店 2430 ((-66390.74 26804.54, -66391.2…
6 セブンイレブン佐賀駅前中央2丁目店 2422 ((-64759.35 29882.63, -64759.8…
7 ミニストップ佐賀田代2丁目店 2417 ((-63594.61 27884.37, -63595.0…
8 セブンイレブン佐賀高木瀬東2丁目店 2407 ((-65423.24 31483.4, -65423.72…
9 セブンイレブン佐賀南本庄店 2364 ((-66338.53 26854.67, -66339.0…
10 セブンイレブン佐賀医大通り店 2350 ((-67705.79 31232.43, -67706.2…
# ℹ 100 more rows
この結果を使って、商圏人口(population)で色をつけて地図にしてみましょう。
ggplot() +
geom_sf(data = saga_city) +
geom_sf(aes(fill = population), data = market) +
scale_fill_distiller(direction = 1) +
theme_minimal()
地図を見ると、人口密度の高い市中心部の店舗の商圏人口が多い傾向にあることがわかります。
佐賀市内のコンビニについて、店舗ごとの商圏人口の分布はどのようになっているでしょうか? ggplot
を使って、ヒストグラムを書いてみましょう。
ggplot(aes(x = population), data = market) +
geom_histogram(boundary = 0, binwidth = 250) +
theme_minimal()
Warning: Removed 2 rows containing non-finite outside the scale range
(`stat_bin()`).
佐賀市のコンビニ立地においては、半径350 mの範囲に1,000人の人口がいるかどうかで、店舗立地の傾向が分かれそうな気がします。 ということで、商圏人口1,000人を閾値にして、地図を書いてみましょう。 mutate
関数を使って、人口が1,000人よりも多ければTRUE
、そうでなければFALSE
となるようなデータ列locationを作成し、locationを使って色を塗り分けた地図を描いてみます。
<- market |>
market mutate(location = population > 1000)
ggplot() +
geom_sf(data = saga_city) +
geom_sf(aes(fill = location), data = market) +
theme_minimal()
先ほどの地図よりも中心部と周辺部の店舗による商圏人口の違いがはっきりとわかる地図になりました。 周辺部のコンビニ店舗は、店舗周辺の人口は少ないけれども、幹線道路沿いであるなど、別の立地要因が影響している可能性があります(道路交通量の分析は今後の課題ということにしておきます)。
新規出店場所の検討
次に、新しく佐賀市内のどこにコンビニを出店可能かについて考えてみましょう。 ここでは、地域メッシュ統計のデータを使ってみることにします。
まず、佐賀市の250 mメッシュのポリゴンデータを作成します。 いろいろな方法がありますが、ここではjpgrid
パッケージのgeometry_to_grid
関数とgrid_as_sf
関数を使った方法を紹介します。
この関数jpgrid::geometry_to_grid
を使えば、行政境界などのsfオブジェクトを、地域メッシュコードに変換することができます。 また、jpgrid::grid_as_sf
により、地域メッシュコード(gridクラス)を含むデータを、sfオブジェクトに変換することが可能です。
佐賀市のポリゴンデータを使って、佐賀市の市域にかかる4次メッシュ(250 mメッシュ)からなるポリゴンデータを作成します。 jpgrid::geometry_to_grid
関数は長さ1のリストを返すので、dplyr::first
関数を使ってリストの要素を取り出しています。 また、後で統計データを結合するためのmeshcodeを文字列型のデータとして作っておきます。
<- saga_city |> st_transform(6668) |>
saga_grid geometry_to_grid(grid_size = '250m') |>
first() |> grid_as_sf(crs = 6668) |>
st_transform(6670) |>
mutate(meshcode = as.character(grid))
できたメッシュデータを佐賀市のポリゴンに重ねて地図にして、確認してみましょう。
ggplot() +
geom_sf(data = saga_city) +
geom_sf(data = saga_grid) +
theme_minimal()
さて、出店候補地を調べるわけですが、ここでは考える出店地域を市南部の平野に限定しましょう。 たまたま、1次メッシュ区画の境界が平地の山地の境目付近にありますので、ここで区切れば良さそうです。 ということで、1次メッシュ区画が4930であるメッシュだけを抽出しましょう(1次メッシュごとに統計データが公開されているので、データダウンロードが1回で済むというメリットもあります)。
ちなみに、九州周辺の1時メッシュコードは下図のようになっています。 佐賀県と重なっているのは、4929、4930、5029、5030の4つで、そのうち佐賀市にかかるのは、 4930と5030の2つです。
5次メッシュコードは11桁の数字ですが、その先頭4桁は対応する1次メッシュコードになっていますので、 メッシュコードの先頭が4930から始まっているメッシュだけをフィルタリングします。 stringr::str_starts
関数は文字列の先頭が第1引数の文字列が、第2引数のパターンから始まる場合にTRUE
を、そうでない場合にFALSE
を返します。
<- saga_grid |>
saga_grid filter(str_starts(meshcode, "4930"))
もう一度地図を表示して、フィルタリングがうまくいったか確認しましょう。
ggplot() +
geom_sf(data = saga_city) +
geom_sf(data = saga_grid) +
theme_minimal()
次に地域メッシュ統計の人口データをダウンロードします。 e-Statから、1次メッシュコード4930の人口データを取得しましょう。
> 統計データダウンロード - [国勢調査]→[2020年]→[5次メッシュ(250mメッシュ)]→[人口及び世帯]→[M4930]→CSVボタンをクリックしダウンロード
「都道府県で絞込みはコチラ」をクリックし、「41 佐賀県」の左のチェックボックスにチェック☑️をいれ「選択」ボタンをクリックすると、佐賀県にかかる4つの1次メッシュだけが表示されるので、そこから「M4930」を選ぶと良いと思います。
tblT001142Q4930.zipがダウンロードされたら、これを展開して、できたファイルtblT001142Q4930.txtだけを、プロジェクトのdataフォルダに移動してください。
さて、これをreadr::read_csv
関数で読み込むわけですが、例によって、2回に分けて読み込みます。 以下のコードでは、1回目はヘッダ情報のみ、2回目はデータ部分のみ読み込んでいます。
<-
col_names read_csv('data/tblT001142Q4930.txt', n_max = 1,
col_names = FALSE, col_types = 'c') |>
as.character()
<-
mesh read_csv('data/tblT001142Q4930.txt', skip = 2,
col_names = col_names, na = c("", "*"))
Rows: 32215 Columns: 54
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): GASSAN
dbl (53): KEY_CODE, HTKSYORI, HTKSAKI, T001142001, T001142002, T001142003, T...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
今回使用するのは、メッシュ総人口のデータだけですので、dplyr::select
関数で抜き出しましょう。 同時に、人口データの列名をpopulationに、またmeshcodeを数値型から文字列型に変更しています。
<- mesh |>
mesh select(meshcode = KEY_CODE,
population = T001142001) |>
mutate(meshcode = as.character(meshcode))
glimpse(mesh)
Rows: 32,215
Columns: 2
$ meshcode <chr> "4930010022", "4930010023", "4930010024", "4930010034", "49…
$ population <dbl> 3, 57, 70, 3, 24, 92, 2, 22, 37, 5, 1, 86, 60, 18, 61, 1, 3…
このmeshをsaga_gridに結合します。dplyr::left_join
を使います。
<-
saga_grid |> left_join(mesh, by = join_by(meshcode))
saga_grid glimpse(saga_grid)
Rows: 3,319
Columns: 4
$ grid <grd250> 4930526634, 4930526643, 4930526644, 4930526733, 49305267…
$ geometry <POLYGON [m]> POLYGON ((-62684.14 15679.9..., POLYGON ((-62392.58…
$ meshcode <chr> "4930526634", "4930526643", "4930526644", "4930526733", "49…
$ population <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
メッシュ人口の地図を描いてみましょう。 メッシュ統計の利点の1つは、各メッシュの大きさがほぼ同じなので、人口密度を計算しなくても、それぞれのメッシュの人口の大小を直接比較できることです。
ggplot() +
geom_sf(aes(fill = population), data = saga_grid, color = NA) +
scale_fill_distiller(direction = 1) +
theme_minimal()
さて、佐賀市におけるコンビニの成立条件は、半径350 mの範囲に1,000の人口がいることでした。 これをメッシュ人口に置き換えると、メッシュの面積がおよそ2502 m2、半径350 mの円の面積がπ× 3502 m2ですから、メッシュ区画が3×3=9個並んだ範囲において、
1000 * 250^2 * 9 / (pi * 350^2 )
[1] 1461.627
およそ1,500人の人口がいることに相当します。 それでは、9区画の人口が1,500人を超えるメッシュ区画を抽出してみましょう。
Rでも他のプログラミング言語と同様に、自前の関数を定義できます。 その基本的な書式は、以下のとおりです。
<- function(variables){
new_function # ここに処理を書く
}
ここでは、例としてベクトルの平均を計算する関数average
を自作してみましょう。
<- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
x <- function(x){
average <- length(x)
n <- sum(x)
s / n
s
}average(x)
[1] 5.5
以下では、周辺9メッシュ合計の人口を計算する関数compute_neighbor_pop
を定義しています。 そして、その関数をpurrr::map_int
関数によって全てのメッシュについて繰り返し計算しています。 計算した結果は、neighbor_pop列に格納しています。
<- function(grid){
compute_neighbor_pop <- grid |>
neighbor_grid grid_neighborhood(n = 0:1, type = 'moore') |>
first()
|>
saga_grid filter(grid %in% neighbor_grid) |>
pull(population) |>
sum(na.rm = TRUE)
}<- saga_grid |>
saga_grid mutate(neighbor_pop = map_int(grid, compute_neighbor_pop))
次に、mutate
関数を使って、人口が1,000人を超えるかどうかをTRUEかFALSEかで格納したmarketというデータ列を作成します。
<- saga_grid |>
saga_grid mutate(candidate = neighbor_pop > 1500)
candidateでメッシュの色を塗り分けた地図を表示してみましょう。 その地図に、現在のコンビニの立地を重ねます。
ggplot() +
geom_sf(aes(fill = candidate), data = saga_grid, color = NA) +
geom_sf(data = buffer) + ylim(16000, 37000) +
theme_minimal()
人口が1,000人を超えていて、なおかつ既存のコンビニが出店していないメッシュが、まずは新規出店の候補地となるのではないでしょうか。 しかし、実際の出展計画を考える際にはこれだけの分析では不十分だと考えられます。 まずは、実際にその付近の土地や建物の利用可能性を確認する必要があります。 また、夜間人口だけでなく、昼間の周辺人口(職場や学校など)、敷地前面道路の交通量など、周辺の土地利用や建築に関わる規制など、多くの要素を考慮に入れる必要があるでしょう。
実際にローソンの出店ガイドラインをみると、「第一種低層住居専用地域、工業専用地域には原則出店できません」との記述があります。 そこで最後に、佐賀市中心部の用途地域ついて確認してみましょう。
都市における住居、商業、工業といった土地利用は、似たようなものが集まっていると、それぞれにあった環境が守られ、効率的な活動を行うことができます。しかし、種類の異なる土地利用が混じっていると、互いの生活環境や業務の利便が悪くなります。 そこで、都市計画では都市を住宅地、商業地、工業地などいくつかの種類に区分し、これを「用途地域」として定めています。 それぞれの用途地域に応じて、建てられる建築物の用途や規模が決められています。
佐賀市の最新の用途地域については、佐賀市ウェブサイトおよびそこからリンクされているぐるっとさがナビから確認することができます。 また、用途地域のポリゴンデータが国土数値情報から入手可能です。
こちらから、2019年(令和元年)の佐賀県の用途地域データ(A29-19_41_GML.zip)をダウンロードしてください。 ダウンロードフォルダにあるzipファイルを展開して、できたフォルダ(A29-19_41_GML)を、プロジェクトのdataフォルダに移動してください。
シェープファイルをsf::st_read
関数で読み込みます。 日本語の文字化けを解消するために、文字コードをoptions引数で指定しています。
<-
zoning st_read('data/A29-19_41/01-03_GeoJSON形式/A29-19_41201.geojson') |>
filter(A29_005 %in% c("第一種低層住居専用地域", "工業専用地域")) |>
st_transform(6670)
Reading layer `A29-19_41201' from data source
`/Users/kzktmr/Documents/LectureNotes/data/A29-19_41/01-03_GeoJSON形式/A29-19_41201.geojson'
using driver `GeoJSON'
Simple feature collection with 188 features and 8 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 130.2618 ymin: 33.20996 xmax: 130.3746 ymax: 33.32395
Geodetic CRS: WGS 84
佐賀市の用途地域データから、用途地域が第一種低層住居専用地域、工業専用地域であるポリゴンだけを抜き出すために、filter
関数でフィルタリングを行っています。 そして、平面直角座標系II系に投影変換しています。
先ほどの地図に、用途地域(第一種低層住居専用地域、工業専用地域)のポリゴンを重ねて表示してみましょう。
ggplot() +
geom_sf(aes(fill = candidate), data = saga_grid, color = NA) +
geom_sf(aes(fill = A29_005), data = zoning) +
geom_sf(data = buffer) + ylim(16000, 37000) +
theme_minimal()
なんとなく、どのあたりに出店したら良さそうか、見えてきたでしょうか。