14 ジオマーケティング

必須パッケージ

  • この章では、以下のパッケージが必要である(revgeoもインストールされている必要がある)。
  • 必要なデータ、順次ダウンロード予定

読者の利便性と再現性を確保するため、ダウンロードしたデータを spDataLarge パッケージで公開している。

14.1 イントロダクション

この章では、パート I とパート II で学んだスキルを特定のドメインに適用する方法を示す。ジオマーケティング (ロケーション分析 やロケーションインテリジェンスとも呼ばれることがある)である。 研究・実用化されている分野は幅広い。 その典型的な例が、新しい店舗をどこに置くかということである。 ここでの目的は、最も多くの訪問者を集め、最終的に最も多くの利益を上げることである。 また、例えば新しい医療サービスをどこに配置するかなど、この技術を公共の利益のために利用できる非商業的なアプリケーションも多い (Tomintz, Clarke, and Rigby 2008)

ロケーション分析の基本は人である。 特に時間やその他のリソースを費やす可能性が高い場所である。 興味深いことに、エコロジーの概念やモデルは、店舗立地分析に使われるものと非常によく似ている。 動物や植物は、空間的に変化する変数に基づいて、特定の「最適な」場所でそのニーズを最もよく満たすことができる(Muenchow et al. (2018); Chapter 15 も参照)。 これはジオコンピューティングやGISサイエンス全般の大きな強みである。 コンセプトや手法は他の分野にも転用可能である。 例えば、ホッキョクグマは気温が低く、餌(アザラシやアシカ)が豊富な北緯を好む。 同様に、人間は特定の場所に集まる傾向があり、北極の生態学的ニッチに類似した経済的ニッチ(そして高い地価)を作り出する。 ロケーション分析の主な作業は、利用可能なデータに基づいて、特定のサービスにとってそのような「最適な場所」がどこであるかを見つけ出すことである。 典型的なリサーチクエスチョンとしては、以下のようなものがある。

  • ターゲット層はどこに住んでいて、どのエリアによく行くのか?
  • 競合する店舗やサービスはどこにあるのか?
  • 特定の店舗にどれくらいの人が行きやすいか?
  • 既存のサービスは、市場の潜在力を過大に、あるいは過小に開拓していないか?
  • ある企業の特定地域における市場シェアはどのくらいか?

本章では、実データに基づく仮想的なケーススタディに基づき、ジオコンピューティングがそのような疑問に答えることができることを示す。

14.2 ケーススタディ: ドイツの自転車店

あなたがドイツで自転車ショップのチェーンを始めたとする。 店舗は、できるだけ多くの潜在顧客がいる都市部に配置すること。 さらに、仮定の調査(この章のために考案されたもので、商業利用はできない!)によると、独身の若い男性(20歳から40歳)が貴社の製品を購入する可能性が最も高いということである:これがターゲット層である。 あなたは、何店舗も出店できる十分な資金をお持ちの幸運な立場にある。 しかし、どこに配置すればいいのだろうか? コンサルティング会社(ジオマーケティング アナリストを雇っている)は、このような質問に答えるために喜んで高い料金を取るだろう。 幸い、オープンデータ( )やオープンソースソフトウェア( )の力を借りれば、私たち自身でそれを行うことができる。 以下の章では、本書の最初の章で学んだテクニックを、サービスロケーション解析の一般的なステップにどのように応用できるかを紹介する。

  • ドイツの国勢調査の入力データを整頓する ( Section 14.3 )
  • 集計された国勢調査データをラスタに変換 オブジェクト ( Section 14.4 )
  • 人口密度の高い都市圏の特定 ( Section 14.5 )
  • これらの地域の詳細な地理データ(OpenStreetMap , with osmdata )をダウンロードする ( Section ?? ) 。
  • 地図代数を用いて異なる場所の相対的な望ましさをスコアリングするためのラスタ を作成する ( Section 14.7 )

これらのステップは、特定のケーススタディに適用されたが、店舗立地や公共サービス提供の多くのシナリオに一般化することが可能である。

14.3 入力データを整頓

ドイツ政府は、1kmまたは100mの解像度でグリッド化された国勢調査データを提供している。 次のコードは、1kmのデータをダウンロードし、解凍し、読み込むものである。 なお、census_despDataLarge パッケージ ( data("census_de", package = "spDataLarge" ) からも入手可能である。

download.file("https://tinyurl.com/ybtpkwxz", 
              destfile = "census.zip", mode = "wb")
unzip("census.zip") # unzip the files
census_de = readr::read_csv2(list.files(pattern = "Gitter.csv"))

census_de オブジェクトは、ドイツ全土の30万以上のグリッドセルについて、13の変数を含むデータフレームである。 私たちの仕事では、これらのサブセットだけが必要である。東経(x)と北緯(y)、住民数(人口;pop)、平均年齢(mean_age)、女性の割合(women)、平均世帯人員(hh_size)である。 これらの変数は、以下のコードチャンクで選択され、ドイツ語から英語に名前が変更され、 Table 14.1 に要約される。 さらに、mutate_all() は、値 -1 と -9 (不明を意味する) を NA に変換するために使用される。

# pop = population, hh_size = household size
input = dplyr::select(census_de, x = x_mp_1km, y = y_mp_1km, pop = Einwohner,
                      women = Frauen_A, mean_age = Alter_D,
                      hh_size = HHGroesse_D)
# set -1 and -9 to NA
input_tidy = mutate_all(input, list(~ifelse(. %in% c(-1, -9), NA, .)))
TABLE 14.1: Categories for each variable in census data from Datensatzbeschreibung…xlsx located in the downloaded file census.zip (see Figure 13.1 for their spatial distribution).
class Population % female Mean age Household size
1 3-250 0-40 0-40 1-2
2 250-500 40-47 40-42 2-2.5
3 500-2000 47-53 42-44 2.5-3
4 2000-4000 53-60 44-47 3-3.5
5 4000-8000 >60 >47 >3.5
6 >8000

14.4 国勢調査ラスタを作成

前処理を行った後、ラスタスタックまたはブリック に変換することができる(Section 2.3.43.3.1 を参照)。 rasterFromXYZ() を使えば簡単である。 入力データフレームが必要で、最初の2列は規則的なグリッド上の座標を表している。 残りの列(ここでは、pop , women , mean_age , hh_size )は、ラスタブリックレイヤ( Figure 14.1 ; github リポジトリの code/14-location-jm.R も参照)の入力として機能する。

input_ras = rasterFromXYZ(input_tidy, crs = st_crs(3035)$proj4string)
input_ras
#> class : RasterBrick
#> dimensions : 868, 642, 557256, 4 (nrow, ncol, ncell, nlayers)
#> resolution : 1000, 1000 (x, y)
#> extent : 4031000, 4673000, 2684000, 3552000 (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=52 +lon_0=10
#> names       :  pop, women, mean_age, hh_size 
#> min values  :    1,     1,        1,       1 
#> max values  :    6,     5,        5,       5
Note that we are using an equal-area projection (EPSG:3035; Lambert Equal Area Europe), i.e., a projected CRS where each grid cell has the same area, here 1000 x 1000 square meters. Since we are using mainly densities such as the number of inhabitants or the portion of women per grid cell, it is of utmost importance that the area of each grid cell is the same to avoid ‘comparing apples and oranges.’ Be careful with geographic CRS where grid cell areas constantly decrease in poleward directions (see also Section 2.4 and Chapter 7).
Gridded German census data of 2011 (see Table 13.1 for a description of the classes).

FIGURE 14.1: Gridded German census data of 2011 (see Table 13.1 for a description of the classes).

次に、input_ras に格納されているラスタの値を、Section 14.2 で述べた調査に従って、Section 4.3.3 で紹介した raster 関数 reclassify() を用いて再分類している。 母集団データの場合、クラスの平均値を用いて数値データ型に変換する。 ラスタセルは、値1(「クラス1」のセルが3~250人の住民を含む)の場合は127人、値2(250~500人の住民を含む)の場合は375人、と仮定される( Table 14.1 を参照)。 これらのセルには8000人以上の人が含まれているため、「クラス6」のセル値には8000人の住民が選ばれた。 もちろん、これは真の母集団の近似値であり、正確な値ではない。93 しかし、大都市圏を定義するには十分なレベルである(次章参照)。

総人口の絶対推計を表す変数 pop とは対照的に、残りの変数は、調査で使用されたウェイトに対応するウェイトに分類し直した。 例えば、変数 women のクラス1は、人口の0~40%が女性である地域を表す。 は、ターゲット層が男性であるため、比較的高いウェイトである3に分類し直した。 同様に、若年層や単身世帯の割合が高い層は、高いウェイトを持つように分類し直した。

rcl_pop = matrix(c(1, 1, 127, 2, 2, 375, 3, 3, 1250, 
                   4, 4, 3000, 5, 5, 6000, 6, 6, 8000), 
                 ncol = 3, byrow = TRUE)
rcl_women = matrix(c(1, 1, 3, 2, 2, 2, 3, 3, 1, 4, 5, 0), 
                   ncol = 3, byrow = TRUE)
rcl_age = matrix(c(1, 1, 3, 2, 2, 0, 3, 5, 0),
                 ncol = 3, byrow = TRUE)
rcl_hh = rcl_women
rcl = list(rcl_pop, rcl_women, rcl_age, rcl_hh)

なお、リスト中の再分類行列の順序は、input_ras の要素と同じになるようにした。 例えば、最初の要素はどちらの場合も母集団に対応する。 その後、for -loop 、再分類行列を対応するラスタレイヤーに適用する。 最後に、以下のコードで、reclass のレイヤが input_ras のレイヤと同じ名前であることを確認する。

reclass = input_ras
for (i in seq_len(nlayers(reclass))) {
  reclass[[i]] = reclassify(x = reclass[[i]], rcl = rcl[[i]], right = NA)
}
names(reclass) = names(input_ras)
reclass
#> ... (full output not shown)
#> names       :  pop, women, mean_age, hh_size 
#> min values  :  127,     0,        0,       0 
#> max values  : 8000,     3,        3,       3

14.5 大都市圏を定義

大都市圏とは、50万人以上が住む20km2のピクセルと定義している。 この粗い解像度のピクセルは、Section 5.3.3 で紹介したように、aggregate() 、速やかに作成することができる。 以下のコマンドは、引数 fact = 20、結果の解像度を20倍にしている(元のラスタの解像度が1km2であったことを思い出してみよう)。

pop_agg = aggregate(reclass$pop, fact = 20, fun = sum)

次のステージは、50万人以上のセルだけを残すことである。

summary(pop_agg)
#>             pop
#> Min.        127
#> 1st Qu.   39886
#> Median    66008
#> 3rd Qu.  105696
#> Max.    1204870
#> NA's        447
pop_agg = pop_agg[pop_agg > 500000, drop = FALSE] 

これをプロットすると、8つの大都市圏(Figure 14.2)が見えてくる。 各領域は、1つ以上のラスタセルで構成される。 1つのリージョンに属するすべてのセルを結合できればコマンドは、 ラスタclump() である。 その後、rasterToPolygons() でラスタオブジェクトを空間ポリゴンに変換し、st_as_sf()sf -オブジェクトに変換する。

polys = pop_agg %>% 
  clump() %>%
  rasterToPolygons() %>%
  st_as_sf()

polysclumps 列は、各ポリゴンがどの大都市圏に属しているかを示し、 ポリゴンを首尾一貫した単一の地域に分解(ディゾルブ)するために使用される(Section 5.2.7 も参照)。

metros = polys %>%
  group_by(clumps) %>%
  summarize()

他の列を入力として与えなければ、summarize() はジオメトリを分解(ディゾルブ)するだけである。

The aggregated population raster (resolution: 20 km) with the identified metropolitan areas (golden polygons) and the corresponding names.

FIGURE 14.2: The aggregated population raster (resolution: 20 km) with the identified metropolitan areas (golden polygons) and the corresponding names.

その結果、自転車店に適した8つの都市圏( Figure 14.2 ;図の作成については code/14-location-jm.R も参照)が、まだ名前が出てこない。 リバースジオコーディング のアプローチでこの問題を解決することができる。 座標が与えられると、逆ジオコーディングにより対応する住所が求められる。 その結果、各都市圏の重心 座標を抽出することで、リバースジオコーディング API の入力とすることができる。 revgeo パッケージは、OpenStreetMap、Google Maps、Bing 用のオープンソース Photon geocoder にアクセスするためのものである。 デフォルトでは、Photon APIを使用する。 revgeo::revgeo() は地理座標(緯度/経度)しか受け付けない。したがって、最初の要件は、大都市ポリゴンを適切な座標参照系 (緯度経度) に持っていくことである。 Chapter 7

metros_wgs = st_transform(metros, 4326)
coords = st_centroid(metros_wgs) %>%
  st_coordinates() %>%
  round(4)

framerevgeocode()output オプションとして選択すると、通り名、家屋番号、都市名など、場所を示すいくつかの列を持つ data.frame が返される。

library(revgeo)
metro_names = revgeo(longitude = coords[, 1], latitude = coords[, 2], 
                     output = "frame")

読者が全く同じ結果を使用できるようにするため、metro_names オブジェクトとして spDataLarge に入れている。

TABLE 14.2: Result of the reverse geocoding.
city state
Hamburg Hamburg
Berlin Berlin
Wülfrath North Rhine-Westphalia
Leipzig Saxony
Frankfurt am Main Hesse
Nuremberg Bavaria
Stuttgart Baden-Württemberg
Munich Bavaria

全体として、私たちは city 列が大都市名 ( Table 14.2 ) として機能していることに満足している。例外は、Wülfrath が Düsseldorf の大領域に属していることで ある。 したがって、Wülfrath を Düsseldorf ( Figure 14.2 ) に置き換える。 ü のようなウムラウトは、例えば opq() を使って大都市圏のバウンディングボックスを決定する場合(後述)、後々トラブルになる可能性があるため、避けているのである。

metro_names = dplyr::pull(metro_names, city) %>% 
  as.character() %>% 
  ifelse(. == "Wülfrath", "Duesseldorf", .)

14.6 注目点

osmdata パッケージは、OSM データへの使いやすいアクセスを提供する ( Section 8.2 も参照)。 ドイツ全土の店舗をダウンロードするのではなく、定義された大都市圏にクエリを限定することで、計算負荷を軽減し、関心のあるエリアのみの店舗位置を提供している。 この後のコードチャンクは、以下のようないくつかの関数を用いてこれを行う。

  • map() ( lapply()tidyverse 相当)。これは、OSM クエリー関数 opq() ( Section 8.2 参照) のバウンディングボックス を定義する、8つの大都市名すべてを繰り返し処理するものである。
  • add_osm_feature() で、キー値が shop の OSM 要素を指定する (共通のキー:値のペアの一覧は wiki.openstreetmap.org を参照)。
  • osmdata_sf()、これは OSM データを空間オブジェクト(クラス sf )に変換するものである。
  • while() 、1回目のダウンロードに失敗すると繰り返し(今回は3回)ダウンロードを試みる94 このコードを実行する前に: 約2GBのデータをダウンロードすることを考慮してみよう。 時間とリソースを節約するために、shops という名前の出力を spDataLarge に入れることとした。 お使いの環境で利用できるようにするには、spDataLarge パッケージがロードされていることを確認するか、data("shops", package = "spDataLarge") を実行してみよう。
shops = map(metro_names, function(x) {
  message("Downloading shops of: ", x, "\n")
  # give the server a bit time
  Sys.sleep(sample(seq(5, 10, 0.1), 1))
  query = opq(x) %>%
    add_osm_feature(key = "shop")
  points = osmdata_sf(query)
  # request the same data again if nothing has been downloaded
  iter = 2
  while (nrow(points$osm_points) == 0 & iter > 0) {
    points = osmdata_sf(query)
    iter = iter - 1
  }
  points = st_set_crs(points$osm_points, 4326)
})

定義された大都市圏に店舗がないことはまずありえない。 次の if の条件は、各リージョンに少なくとも1つのショップがあるかどうかをチェックするだけである。 その場合は、該当する地域のショップを再度ダウンロードすることを勧める。

# checking if we have downloaded shops for each metropolitan area
ind = map(shops, nrow) == 0
if (any(ind)) {
  message("There are/is still (a) metropolitan area/s without any features:\n",
          paste(metro_names[ind], collapse = ", "), "\nPlease fix it!")
}

各リスト要素 ( sf データフレーム) が同じ列を持つことを確認するために、osm_idshop の列だけを、別の map ループの助けを借りて保持する。 OSMの投稿者は、データを収集する際に同じように細心の注意を払っているわけではないので、これは当たり前のことではない。 rbind 最後に、すべてのショップを1つの大きな sf オブジェクトにまとめる。

# select only specific columns
shops = map(shops, dplyr::select, osm_id, shop)
# putting all list elements into a single data frame
shops = do.call(rbind, shops)

単純に map_dfr() とした方が簡単だっただろう。 残念ながら、今のところ、sf のオブジェクトと調和して動作しない。 注: shops は、spDataLarge パッケージで提供される。

あとは、空間点オブジェクトをラスタに変換するだけである( Section 6.4 参照)。 sf オブジェクト shops は、reclass オブジェクトと同じパラメータ(寸法、解像度、CRS )を持つラスタ に変換される。 重要なのは、ここで count() 関数を用いて、各セルのショップ数を算出していることだ。

If the shop column were used instead of the osm_id column, we would have retrieved fewer shops per grid cell. This is because the shop column contains NA values, which the count() function omits when rasterizing vector objects.

そのため、後続のコードチャンクの結果は、店舗密度(店舗/km2)の推定値となる。 st_transform() は、両入力の CRS が一致するように、 rasterize() の前に使用される。

shops = st_transform(shops, proj4string(reclass))
# create poi raster
poi = rasterize(x = shops, y = reclass, field = "osm_id", fun = "count")

他のラスタレイヤ(人口、女性、平均年齢、世帯人員)と同様、poi ラスタは4つのクラスに再分類される(Section 14.4 参照)。 クラス間隔の定義は、ある程度恣意的に行われるものである。 均等割、分位割、固定値などを使用することができる。 ここでは、クラス内分散を最小化する Fisher-Jenks 自然休息法を選択し、その結果を再分類行列の入力とする。

# construct reclassification matrix
int = classInt::classIntervals(values(poi), n = 4, style = "fisher")
int = round(int$brks)
rcl_poi = matrix(c(int[1], rep(int[-c(1, length(int))], each = 2), 
                   int[length(int)] + 1), ncol = 2, byrow = TRUE)
rcl_poi = cbind(rcl_poi, 0:3) 
# reclassify
poi = reclassify(poi, rcl = rcl_poi, right = NA) 
names(poi) = "poi"

14.7 適当な場所を特定

すべてのレイヤを結合する前に残っている唯一のステップは、poireclass のラスタスタックに追加し、そこから人口レイヤを削除することである。 後者の理由は2つある。 まず、大都市圏、つまりドイツの他の地域に比べて人口密度が平均的に高い地域はすでに定義されている。 第二に、特定のキャッチメントエリア内に多くの潜在顧客がいることは有利であるが、数が多いだけでは、実際には望ましいターゲットグループを表していない可能性がある。 例えば、タワーマンションは人口密度が高い地域であるが、高価なサイクル部品の購買力が高いとは限らない。 これは、相補関数 addLayer()dropLayer() によって実現される。

# add poi raster
reclass = addLayer(reclass, poi)
# delete population raster
reclass = dropLayer(reclass, "pop")

他のデータサイエンス・プロジェクトと同様、これまでのところ、データの検索と「整理」が全体の作業負荷の多くを占めている。 きれいなデータであれば、最後のステップであるすべてのラスタのレイヤを合計して最終的なスコアを計算することも、1行のコードで実現できる。

# calculate the total score
result = sum(reclass)

例えば、9以上のスコアは、自転車ショップを配置できるラスタセルを示す適切な閾値かもしれない(Figure 14.3 ; code/14-location-jm.R も参照)。

FIGURE 14.3: Suitable areas (i.e., raster cells with a score > 9) in accordance with our hypothetical survey for bike stores in Berlin.

14.8 ディスカッションと次のステップ

今回紹介したアプローチは、GIS の規範的な使い方の典型的な例である (Longley 2015)。 調査データと専門家による知識・仮定(大都市圏の定義、クラス間隔の定義、最終的なスコア閾値の定義)を組み合わせている。 このアプローチは、科学的な研究よりも、他の情報源と比較すべき、自転車店に適した地域のエビデンスに基づく指標を提供する応用分析に適している。 アプローチにいくつかの変更を加えることで、分析結果を改善することができる。

  • 最終的なスコアの算出には均等なウェイトを用いたが、世帯規模など他の要因も、女性の割合や平均年齢と同様に重要である可能性がある。
  • 全てのポイントオブインタレスト を使用したが、DIY、ハードウェア、自転車、釣り、ハンティング、バイク、アウトドア、スポーツショップなど、バイクショップに関連するもののみ(ショップ値の範囲は OSM Wiki で確認可能)、より洗練された結果を得ることができたかもしれない
  • より高い解像度のデータにより、出力が向上する場合がある(演習参照)
  • 私たちは限られた変数のみを使用し、 INSPIRE geoportal やOpenStreetMapのサイクリングロードのデータなど、他の情報源からのデータは分析を豊かにするかもしれない(Section 8.2 も参照のこと)。
  • 男性比率と単身世帯の関係などの相互作用は考慮されていない。

つまり、この分析は多方面に拡張できるのである。 とはいえ、ジオマーケティングの文脈の中で、R で空間データを取得し、扱う方法について、第一印象と理解を深めていただけたと思われる。

最後に、今回の分析は、あくまでも適地探しの第一歩に過ぎないということを指摘しておく必要がある。 これまでの調査により、1km四方でバイクショップの立地が可能なエリアを特定した。 その後の分析のステップを踏むことができる。

  • 特定のキャッチメントエリア内の住民の数に基づいて最適な場所を見つける . 例えば、できるだけ多くの人が自転車で15分以内の移動距離で行けるお店であること(キャッチメントエリア ルーティング )。 そのため、店舗から遠ければ遠いほど、実際に店舗を訪れる可能性が低くなることを考慮する必要がある(距離減衰関数)。
  • また、競合他社を考慮するのも良いアイデアだろう。 つまり、選択した場所の近辺にすでに自転車屋がある場合、可能性のある顧客(または販売可能性)を競合他社に分散させる必要がある (Huff 1963; Wieland 2017) .
  • 例えば、アクセスの良さ、駐車場の有無、通行人の希望頻度、大きな窓があることなど、適切かつ手頃な価格の不動産を探す必要がある。

14.9 演習

  1. raster::rasterFromXYZ() を使って、input_tidy をラスタブリックに変換しなさい。 sp::gridded() 機能の助けを借りて、同じことを実現してみてみよう。

  2. 100mセル分解能の住民情報を含むcsvファイルをダウンロードする(https://www.zensus2011.de/SharedDocs/Downloads/DE/Pressemitteilung/DemografischeGrunddaten/csv_Bevoelkerung_100m_Gitter.zip?__blob=publicationFile&v=3)。 解凍後のファイルサイズは1.23GBとなるので注意。 R に読み込むには、readr::read_csv を使用する。 私のマシン(16GB RAM)では30秒かかる。 data.table::fread() はさらに高速で、 クラスのオブジェクトを返すだろう。 data.table() as.tibble() を使ってティブルに変換する。 住民ラスタを作成し、セル解像度1kmでこれを集計し、クラス平均値を用いて作成した住民ラスタ(inh)との差を比較しなさい。

  3. 私たちの自転車ショップでは、主に高齢者に電動自転車を販売しているとする。 年齢ラスタを適宜変更し、残りの解析を繰り返し、元の結果と比較しなさい。