R でポリゴンに代表値付与（点データ）

地価、傾斜方位、土地利用項目、ブランドごとのコンビニ、疾患発症者などの量的変数またはカテゴリ変数を持つ点データあるいはポリゴンデータに対して、特定の街区やグリッドメッシュなどのポリゴンで集計し、代表値がほしい場合があります。

データの内容や目的に応じ、量的変数であれば平均値や中央値、最頻値。カテゴリ変数であれば最頻値。また集計ポリゴンの重心最近傍に位置する地物といった基準で、全地物の統計量あるいは代表地物の値を採用したいことがあるでしょう。

GIS といえば最近では QGIS が使われることが多いですが、上記のような処理は不可能ではありませんが、結構たいへんです。（細かな処理を積み重ねる必要があります。）

そんな処理を R でやるには。

## サンプルデータ

サンプルデータとして、[国土数値情報の小学校区データ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A27-v2_1.html)の平成28年神奈川県のデータ `A27-16_14_GML.zip` より小学校点データ `A27P-16_14.shp` と小学校区ポリゴンデータ `A27-16_14.shp` を使用します。

また、あまりよい例ではありませんが、集計するためにデータとして数値型の `id` を付加しています。（すでに行政コードはありますが、今回のデモでは同一行政区となり、値にばらつきがないので）

こんなデータです。グリッドは 2km 四方。

![サンプルデータ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/215290/066b1ae9-2d87-6f58-d5cd-8da860a8de1d.png)


## 集計
### 平均

```r
library(tidyverse)
library(sf)
grid.mean <- school.points %>%
    select(id) %>%
    aggregate(grid, FUN=mean) %>%
    rename(mean=id)
```

集計するには `aggregate` を使います。上記例では集計対象である `school.points` を `grid` で空間的にグルーピングし、関数 `mean` の結果を得ます。得られる結果は `grid` に集計した値が付加されたデータとなります。

```r
print(grid.mean)
```

```
Simple feature collection with 25 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -30000 ymin: -60000 xmax: -20000 ymax: -50000
Projected CRS: JGD2011 / Japan Plane Rectangular CS IX
First 10 features:
       mean                       geometry
1  167.6000 POLYGON ((-30000 -60000, -2...
2  288.0000 POLYGON ((-28000 -60000, -2...
3  157.2500 POLYGON ((-26000 -60000, -2...
4  170.4000 POLYGON ((-24000 -60000, -2...
5  165.5000 POLYGON ((-22000 -60000, -2...
6  210.0000 POLYGON ((-30000 -58000, -2...
7  154.6667 POLYGON ((-28000 -58000, -2...
8  326.5000 POLYGON ((-26000 -58000, -2...
9   21.0000 POLYGON ((-24000 -58000, -2...
10 211.6667 POLYGON ((-22000 -58000, -2...
```

上記の例では、不要な列を除外するために `select` を使ったり、意味合いが変化したため `rename` で列名を変えたりしていますが、重要なものではなく別になくても構いません。

### 中央値

同様に `FUN=median` とします。

```r
grid.median <- school.points %>%
    select(id) %>%
    aggregate(grid, FUN=median) %>%
    rename(median=id)
```

### 最大値

同様に `FUN=max` とします。

```r
grid.max <- school.points %>%
    select(id) %>%
    aggregate(grid, FUN=max) %>%
    rename(max=id)
```

`aggregate` は集計対象の `sf` オブジェクトの各列に対し `FUN` を実行するため、たとえば、（生徒数データが整備されているとして）「生徒数が最多の小学校の名称を代表値にする」といったことは `aggregate` ではできません。

最大をとる地物を代表地物とするには `st_join` で空間結合をし、 `group_by` でグルーピングしたのち、 `filter` で最大 `id` を持つ地物のみにしています。

基本方針としてはこんな感じです。 `id == max(id)` に該当する地物が複数発生しうる場合は、なんらかの追加条件を加える必要があるかと思います。

```r
grid.max_feat <- grid %>%
    rename(grid_id=id) %>%
    st_join(school.points) %>%
    group_by(grid_id) %>%
    filter(id==max(id)) %>%
    ungroup()
```


### 個数

```r
grid.count <- school.points %>%
    select(id) %>%
    aggregate(grid, FUN=length) %>%
    rename(count=id)
```

### 最頻値

R には最頻を求める関数はないようです。そもそも最頻を求めるために、階級を設定する必要があり、また最大度数をとる階級が複数存在することも普通にあるため、最頻値を適切に求めることは意外と難しいです。

とりあえず、ここではカテゴリ変数などに適する離散値に対する各値の最頻値を求める `mode_value` という関数と、量的変数に適するカーネル密度推定による最大推定密度を取る値を求める `mode_kde` を定義してみました。

```r
mode_value <- function(x){
    return (x %>% table() %>%
        sort() %>% rev() %>%
        names() %>% .[[1]] %>% as.numeric())
}

mode_kde <- function(x, bw='nrd0'){
    if(length(x) == 0){
        return(NA)
    } else if(length(x) == 1){
        return(x[1])
    } else{
        d <- density(x, bw=bw)
        return(d$x[which.max(d$y)])
    }
}
```

これらの定義した関数を使います。

```r
grid.mode <- school.points %>%
    select(id) %>%
    aggregate(grid, FUN=mode_value) %>%
    rename(mode=id)
```

```r
grid.mode_kde <- school.points %>%
    select(id) %>%
    aggregate(grid, FUN=mode_kde) %>%
    rename(mode_kde=id)
```

### そのほか

定義した関数の利用が可能なので、 $\mathrm{median} \pm 2\mathrm{mad}$ あたりで検知した外れ値を除去し、その後平均を取るといったこととかもできます。

## 地理的条件
### グリッド重心最近傍

集計するグリッドの中心（重心）に最も近い点を代表値とするには

```r
# 中心点座標を持ったグリッドを作成
grid.centroid <- grid %>%
    rename(grid_id=id) %>%
    mutate(grid_centroid=st_centroid(geometry))

# グリッド重心との距離を取り、近傍点を求める
# （グリッドに含まれる点であることが前提）
school.points.closest <- school.points %>%
    st_join(grid.centroid) %>%
    filter(!is.na(grid_id)) %>% # drop_na を使うと sf じゃなくなる
    mutate(dist_centroid=st_distance(geometry, grid_centroid, by_element=T)) %>%
    group_by(grid_id) %>%
    filter(dist_centroid==min(dist_centroid)) %>%
    ungroup() %>%
    select(-grid_centroid, -dist_centroid) %>%
    rename(school_id=id) %>%
    st_drop_geometry()

# グリッドと連結
grid.closest <- grid %>%
    left_join(school.points.closest, by=c('id' = 'grid_id'))
```

こんな感じでしょうか。

次は、ポリゴンに対する集計をやってみたいと思います。

## （参考）サンプルデータ詳細

今回使用したサンプルデータです。

```r
library(tidyverse)
library(sf)

# データ読み込み
school.points <- st_read('./A27-16_14_GML/shape/A27P-16_14.shp',
        options='ENCODING=cp932', crs=st_crs(6668),
        stringsAsFactors=F, as_tibble=T) %>%
    st_transform(st_crs(6677)) %>%
    mutate(id=row_number()) %>%
    select(id, A27_001, A27_002, A27_003, A27_004)

school.area <- st_read('./A27-16_14_GML/shape/A27-16_14.shp',
        options='ENCODING=cp932', crs=st_crs(6668),
        stringsAsFactors=F, as_tibble=T) %>%
    st_transform(st_crs(6677)) %>%
    select() %>% # 小学校点データと属性名を揃えるため、敢えて列削除
    st_join(school.points) %>%
    select(id, A27_001, A27_002, A27_003, A27_004)

# グリッド
env <- matrix(c(-30000, -60000, -20000, -60000, -20000, -50000,
        -30000, -50000, -30000, -60000), ncol=2, byrow=T)
env <- st_polygon(list(env))
grid <- env %>%
    st_make_grid(cellsize=2000, crs=st_crs(6677)) %>%
    st_as_sf() %>%
    mutate(id=row_number()) %>%
    select(id) %>%
    rename(geometry=x)

# 描画
p <- ggplot() +
    geom_sf(data=school.area, fill='ivory', color='black', size=0.1) +
    geom_sf(data=grid, fill=NA, color='navy', size=1) +
    geom_sf(data=school.points, color='red', size=1) +
    coord_sf(crs=st_crs(6677), datum=st_crs(6677),
        xlim=c(-30000, -20000), ylim=c(-60000, -50000)) +
    scale_x_continuous(labels=scales::comma) +
    scale_y_continuous(labels=scales::comma)
p

# 画像保存
ggsave('plot.png', plot=p)
```