Plot.jlで日本地図にデータを表示する

# モチベーション
　JuliaのPlots.jlで日本の統計データをマップグラフを表示したい。「Makie.jl使えばいいじゃん」ともおもうが、StatsPlotsなど使い慣れたツールで実現したい。こういうグラフを表示したい。
![1_tUJIEH56_hKfZeB_aa3LPA.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/e4574952-982d-89e2-ef0c-b8e833d62a72.png)

## 地図データの入手
[国土交通省国土地理院](https://www.gsi.go.jp/kankyochiri/gm_jpn.html) にある全レイヤデータをダウンロード（zipファイル）して解凍する。Wikipediaによると地図データは3つのファイルがペアで存在する。
- .shp —シェープ規格：地形情報の本体。 
- .shx —シェープインデックス規格：地形データの前方検索、後方検索を高速にするための位置インデックス。 
- .dbf —属性規格：各シェープに対する縦表形式の属性情報。dBASE IV形式準拠。

shpデータがあれば地図を表示するだけならできるので、このフォーマットを扱うLibirayであるShapefile.jlがあればなんとかなる。

## Shapefile.jl 導入して地図データを読み込む

```Julia
   using Plots, DataFrames, CSV, Shapefile
   coast_table = Shapefile.Table("geodata-jpn-all_u_2_2\\coastl_jpn.shp")
   table_jpn = Shapefile.Table("geodata-jpn-all_u_2_2\\polbnda_jpn.shp") 
```
Table関数はshpデータをテーブルデータとして読み込む。テーブルデータには地域の名称やIDなどを含んだものである。shapes関数はTable内の複数のrowをひとつのポリゴンデータにまとめる。まとまったポリゴンデータはPlotsに渡せばそのまま表示される。例えば上記のcoast_shapeのデータは次のようにして表示が可能。
```Julia
    coast_shape = Shapefile.shapes(coast_table)
    plot(coast_shape,color=:gray)
```
![coast_jpn.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/45e14ef4-4791-aa61-80f2-a4590c8945d6.png)
## Tableのデータをのぞいてみる
TableデータはDataFrameへ簡単に変換できる。
```Julia
Julia> df_pola = DataFrame(table_jpn);
5×9 DataFrame
 Row │ f_code  coc     nam        laa            pop      ypc    adm_code  salb    soc    
     │ String  String  String     String         Int64    Int64  String    String  String 
─────┼────────────────────────────────────────────────────────────────────────────────────
   1 │ FA001   JPN     Hokkai Do  Sapporo Shi    1930496   2014  01100     UNK     JPN
   2 │ FA001   JPN     Hokkai Do  Hakodate Shi    274485   2014  01202     UNK     JPN
   3 │ FA001   JPN     Hokkai Do  Otaru Shi       127224   2014  01203     UNK     JPN
   4 │ FA001   JPN     Hokkai Do  Asahikawa Shi   349057   2014  01204     UNK     JPN
   5 │ FA001   JPN     Hokkai Do  Muroran Shi      91276   2014  01205     UNK     JPN
```
Tableのメタデータをたよりに、該当するrowデータをTableから抽出してShape（複数ならShapes）に渡すことで、Plots.jlが表示できるようになる。日本国土地理院のメタデータでもっとも重要なのが adm_code である。これは**市区町村コード**と呼ばれるもので政府各種統計の地域区分に使われている。[総務省のページ](https://www.soumu.go.jp/denshijiti/code.html)に市区町村コード一覧があるので参照されたい。
- 市区町村コードは5桁の数字と下1桁のパリティ（誤入力防止）で構成されている。
- 本国土地理院のメタデータは上5桁だけを使用している。
- 数字であるが整数Integerではなく数字の文字列Stringである。

以上の点に注意し市区町村コードで当該エリアのShapeを返す関数を考えることができる。

## 市町村コードからShapeを返す関数
```Julia
function selectshape(table_jpn, tg)
    rs = empty(Shapefile.shapes(table_jpn))
    for row in table_jpn
        if typeof(tg) == Regex &&  match(tg, row.adm_code) !== nothing
            push!(rs,Shapefile.shape(row))
        elseif  typeof(tg) == String && tg ==row.adm_code
            push!(rs,Shapefile.shape(row))
        end
    end
    return rs
end
```
正規表現で検索できるようにしたのは、市区町村コードは上2桁で都道府県を表していて、例えば北海道なら01～、東京なら13～で始まるルールを持っている。正規表現を使えば次のような検索ができる。
```Julia
Julia> hokkaido = selectshape(table_jpn, r"^01d*");　#01~は北海道
Julia> plot(hokkaido,alpha=0.3,color=:green)
```
![hokkaido.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/9b52f07e-f8ef-65a4-be95-a8841c41b374.png)

関東のような広域は次のようにしてできよう。
```Julia
kanto = begin
    kanto_rgx = [r"^08d*",r"^09d*",r"^10d*",r"^11d*",r"^12d*",r"^13d*",r"^14d*"]
    # 08~茨木、09~栃木、10~群馬、11~埼玉、12~千葉、13～東京都, 14~神奈川
    [ selectshape(table_jpn, rgx) for rgx in kanto_rgx]
end
pl = plot(xlims=(138,142),ylims=(33.5,37.5))
for (e,c) in zip(kanto,[:yellow,:red,:green,:blue,:orange,:skyblue,:lightgreen])
    plot!(pl, e,color=c,lw=0)
end
display(pl)
```
![kanto.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/3b4bec90-44d4-f8fb-cd0e-06548e6f1f53.png)


## 応用～実際の統計データを表示する

[内閣府のホームページ](https://www5.cao.go.jp/keizai-shimon/kaigi/special/future/keizai-jinkou_data.html)にある1980年と2010年の人口指標（普通出生率：1000人当たり出生数)のマッピング例を再現してみる。

まずはダウンロードしたデータから不必要なデータをカットしてCSVで読み込む。このとき**市区町村コードはTypeをStringとして指定**しなければならない。
```Julia
Julia> birth_da=CSV.read("市区町村別人口1000人当たり年間出生数.csv",
                DataFrame ; types=Dict(:adm_code => String) )
5×6 DataFrame
 Row │ nam            adm_code  y1980    y1990    y2000    y2010   
     │ String         String    Float64  Float64  Float64  Float64 
─────┼─────────────────────────────────────────────────────────────
   1 │ 北海道 札幌市  011002       15.0     10.5      8.5      7.6
   2 │ 北海道 函館市  012025       13.3      9.1      7.5      6.5
   3 │ 北海道 小樽市  012033       10.8      7.1      6.4      5.4
   4 │ 北海道 旭川市  012041       13.9      9.4      8.3      7.4
   5 │ 北海道 室蘭市  012050       13.7      8.1      7.6      6.7
```
plotする
```Julia
n,m=size(birth_da)
pl_80=plot()
pl_10=plot()
year80 = "y1980"
year10 = "y2010"

for i = 1:n
    streetnam = birth_da[i,:adm_code][1:5]
    damap = selectshape(table_jpn, streetnam)
    if length(damap) > 0
        plot!(pl_80,damap, 
            lw=0, 
            clims=(1,25),
            fill_z=birth_da[i, Symbol(year80)],
            fillcolor=:inferno)
        plot!(pl_20,damap, 
            lw=0, 
            clims=(1,25),
            fill_z=birth_da[i, Symbol(year10)], 
            fillcolor=:inferno)
    end
end

L=@layout [a{0.5w} b{0.5w}]
plot(pl_80, pl_20,layout=L,size=(1200,800),dpi=120)
```
これでマップ比較できるようになった。
![市町村出生率.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/379124/82309d7a-1711-0dc2-f762-23fdf300a090.png)

### 今後の課題
- 沖縄と小笠原をきれいにレイアウトしたい。
- GIS地理情報とマッピングできるようにしたい。
- OpenStreetMapX.jlと連携できるようになりたい。
- 図法による正しい投影についても反映していきたい。