SQLだけで Huff Model による商圏分析

# ハフモデルによる商圏分析をSQLでやってみよう！

## 投稿のモチベーション
PostGISとちょっとした工夫で、いろいろなデータ分析ができることを知ってほしいと思い、今回は高額なソフトウェアやサービスとして売られている場合もある商圏分析についてSQLでチャレンジします。

## ハフモデルとは
経済学者 David Huff博士が考案したGISを用いた経済モデルで、店舗の魅力度と距離から、顧客の来店・購入の確率（吸引率）を求めるものです。これがあれば、出店計画の際に取得顧客から売上を予測することも可能です。

ある場所における、１つの店舗の吸引率は、下記の式より求められます。

```math
吸引率 = \frac{魅力度/距離^\alpha}{\sum(魅力度/距離^\alpha)}
\\
\\ここで, \alphaは距離抵抗の係数
```

 - ESRI様解説サイト https://business-map.esrij.com/glossary/2021/
 - PASCO様解説サイト https://www.pasco.co.jp/recommend/word/word036/


ハフモデルのポイントとしては、距離によって魅力度が減衰することです。 
　ディズニーランドとUSJの魅力度が同じとした場合、関西のほとんど人はUSJに行き、関東の人はディズニーランド、東海地方の人はディズニーランドとUSJの半々になるようなイメージです。

# データセット

- 対象とするエリアと店舗
    - 長野県長野市のショッピングセンター （ １２店舗 ）
    <font color="#aaaaaa">台風１９号の水害からの復興を心より願っております。</font>
- テーブル
    - 250m（5次）メッシュ形状とメッシュ人口・世帯数（国勢調査）
        - テーブル名： mesh_l5
        - 列名・型
            - geom geometry(MultiPolygon,4326) -- 5次メッシュの形状
            - key_code bigint -- メッシュコード
            - population_all integer -- メッシュ内人口
            - population_m integer -- メッシュ内男性人口
            - population_f integer -- メッシュ内女性人口
            - households_num integer -- メッシュ内世帯数
    - ショッピングセンターの名称・位置情報・店舗面積 （ 店舗のリストから独自調査 ）
        - テーブル名： shop_list
        - 列名・型
            -  shop_id text -- 店舗id
            -  shop_name text -- 店舗名称
            -  geom geography -- 店舗座標
            -  shop_area integer -- 店舗面積（魅力度）

##### １２店舗の分布状況
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/53918/3c8f2a91-4c1d-0f47-1fc7-7b624e4521e7.png)
© [OpenStreetMap](https://www.openstreetmap.org/copyright) contributors
長野の駅周辺に都市型の店舗が密集しています。


#### ５次メッシュの人口分布
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/53918/4d38836b-ec7a-87f8-fb07-81b79c92ad8e.png)
© [OpenStreetMap](https://www.openstreetmap.org/copyright) contributors


# クエリ

距離抵抗係数(alpha)を 0.5 ~ 2.0 まで0.5刻みでハフモデルによる吸引率を計算します。
 
```sql
CREATE TABLE nagano_trade_area_analysis AS 

WITH 
-- それぞれの店舗と５次メッシュの距離を計算（計算量を減らすため店舗周辺２０kmに限定）
step1 AS (
	SELECT 
	 shop_id
	 , mesh_l5.key_code
	 , ST_Distance(shop_list.geom
		  , ST_Centroid(mesh_l5.geom)::geography )AS distance -- ST_Centroid でメッシュの重心を算出
	 , shop_area
	FROM shop_list, mesh_l5
	WHERE 
		ST_Dwithin(  mesh_l5.geom::geography , shop_list.geom , 20000)
)

SELECT 
  step1.key_code
, shop_id
, alpha
, (shop_area /distance ^ alpha) 
   / SUM(shop_area 
         / (distance ^ alpha) )OVER(PARTITION BY step1.key_code,alpha)
   AS suction_rate
FROM
	 step1 , generate_series(0.5,2.0,0.5) AS alpha ;
```

## 解説

- STEP1では店舗とメッシュの距離の計算を実行しています
  - メッシュ重心の計算 ```ST_Centroid(mesh_l5.geom)::geography )``` 
  - 距離計算 ```ST_Distance(店舗の座標, メッシュ重心の座標)``` 

- ２つめのSELECT句にて、吸引率を一気に算出しています
  - 魅力度の分子 ```(shop_area /distance ^ alpha) ```
  - 魅力度の分母 ```SUM(shop_area / (distance ^ alpha) )OVER(PARTITION BY step1.key_code,alpha)```
  - 距離抵抗係数 0.5 〜 2.0 を０．５毎に発生 ```  generate_series(0.5,2.0,0.5) ```
  
ポイントは、SUM関数をWindow関数として用いていることです。 ```OVER(PARTITION BY step1.key_code,alpha)``` とすることで、メッシュと距離抵抗係数毎のすべての店舗の総和を求めています。 このWindow関数を用いないと、もう１回サブクエリを記述する必要があります。

## 商圏分析
データを抽出し可視化してみます。

### 可視化用データ抽出クエリ
```sql
-- 例） 距離抵抗係数 2.0 のメッシュ毎の吸引率が最大となる値と店舗の抽出
WITH t1 AS ( 
SELECT
      nagano_trade_area_analysis.*    
      , geom
      , rank()OVER( PARTITION BY  mesh_l5.key_code ORDER BY suction_rate DESC )  
FROM nagano_trade_area_analysis  JOIN  mesh_l5 
 ON( nagano_trade_area_analysis.key_code  = mesh_l5.key_code ) 
WHERE alpha = 2.0 ) 
SELECT * FROM t1 WHERE rank = 1;

-- 例） 距離抵抗係数 0.5 のメッシュ毎の吸引率が最大となる値と店舗の抽出
WITH t1 AS ( 
SELECT
      nagano_trade_area_analysis.*    
      , geom
      , rank()OVER( PARTITION BY  mesh_l5.key_code ORDER BY suction_rate DESC )  
FROM nagano_trade_area_analysis  JOIN  mesh_l5 
 ON( nagano_trade_area_analysis.key_code  = mesh_l5.key_code ) 
WHERE alpha = 0.5 ) 
SELECT * FROM t1 WHERE rank = 1;
```

#### 可視化
QGISで表示してみます qgis2threejsプラグインで観察してみます

![dist_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/53918/4892740f-eef7-5488-e449-18316903aeb3.png)

距離抵抗係数が小さいほど広範囲に商圏が広がるので、１つのメッシュの顧客を奪い合う構造になり、結果として各店舗が強く獲得できるエリアは狭くなる傾向が見て取れます。

### 店舗毎の顧客獲得数の計算

今回は距離抵抗係数を２．０の結果を用いて、店舗の毎の獲得人数・世帯数を集計します。

```sql:

SELECT 
 shop_id, alpha 
 , COUNT(nagano_trade_area_analysis.key_code)
 , SUM( population_all*suction_rate ) AS customer_population_all
 , SUM( population_m*suction_rate ) AS customer_population_male 
 , SUM( population_f*suction_rate ) AS customer_population_female
 , SUM( households_num*suction_rate ) AS households_num
FROM 
 nagano_trade_area_analysis 
JOIN  population_mesh_l5 
ON  nagano_trade_area_analysis.key_code = population_mesh_l5.key_code
WHERE alpha = 2.0
GROUP BY 1,2
```
計算結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/53918/c17b7f39-721e-1140-065b-ee4a1eb1684d.png)

長野市では店舗ID１１が強い商圏を持っているという結果になりました。

# 感想
今回はSQLでハフモデルを用いた商圏分析にチャレンジしました。
メッシュ形状の空間データが用意されていれば、GISの処理としては距離計算のみと非常に簡単なものでした。

商圏分析の結果については、地元民のわたしの実感からは離れている感じでした。店舗の面積を魅力度とせってしたところが原因なのかもしれませんし、長野の街の唯一のデーパートである東急が含まれていないことも考えられます。

### 出典

 - 地図でみる統計ダウンロード 
    - 長野市周辺の250ｍメッシュ人口（ [国勢調査 人口等基本集計に関する事項 長野県](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=1&toukeiCode=00200521&toukeiYear=2015&aggregateUnit=Q&serveyId=Q002005112015&statsId=T000876&prefCode=20) ）
    - 長野市周辺の250ｍメッシュ境界([地図で見る統計(統計GIS)  境界データダウンロード](https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=Q&prefCode=20))
 - ショッピングセンター
    -  [日本ショッピングセンター協会 全国都道府県別SC一覧（平成29年）]( http://www.jcsc.or.jp/data/pdf/2018list_of_shopping_centers.pdf )