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

モチベーション

　JuliaのPlots.jlで日本の統計データをマップグラフを表示したい。「Makie.jl使えばいいじゃん」ともおもうが、StatsPlotsなど使い慣れたツールで実現したい。こういうグラフを表示したい。

地図データの入手

国土交通省国土地理院 にある全レイヤデータをダウンロード（zipファイル）して解凍する。Wikipediaによると地図データは3つのファイルがペアで存在する。

• .shp —シェープ規格：地形情報の本体。
• .shx —シェープインデックス規格：地形データの前方検索、後方検索を高速にするための位置インデックス。
• .dbf —属性規格：各シェープに対する縦表形式の属性情報。dBASE IV形式準拠。

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

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

   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のデータは次のようにして表示が可能。

    coast_shape = Shapefile.shapes(coast_table)
    plot(coast_shape,color=:gray)

Tableのデータをのぞいてみる

TableデータはDataFrameへ簡単に変換できる。

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 である。これは市区町村コードと呼ばれるもので政府各種統計の地域区分に使われている。総務省のページに市区町村コード一覧があるので参照されたい。

• 市区町村コードは5桁の数字と下1桁のパリティ（誤入力防止）で構成されている。
• 本国土地理院のメタデータは上5桁だけを使用している。
• 数字であるが整数Integerではなく数字の文字列Stringである。

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

市町村コードからShapeを返す関数

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> hokkaido = selectshape(table_jpn, r"^01d*");　#01~は北海道
Julia> plot(hokkaido,alpha=0.3,color=:green)

関東のような広域は次のようにしてできよう。

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)

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

内閣府のホームページにある1980年と2010年の人口指標（普通出生率：1000人当たり出生数)のマッピング例を再現してみる。

まずはダウンロードしたデータから不必要なデータをカットしてCSVで読み込む。このとき市区町村コードはTypeをStringとして指定しなければならない。

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する

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)

これでマップ比較できるようになった。

今後の課題

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