駅別乗降客数をGoogleEarthで視覚化する

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3116310/bae12e03-9aa8-412c-a65c-d4eef348fd0d.png)
数字がうじゃうじゃしている．こういうのを作ります．途中，色々と調べながら作ったので備忘録として残しておきます．

## 1. データ収集
国土交通省の国土数値情報ダウンロードサイトというところで，駅別の乗降客数のデータが公開されています．

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-S12-2022.html

これをダウンロード・整理して，CSV形式にまとめようと思います．

解凍すると，.shp，.dbf，.prj，.shp，.shxという拡張子のファイルが現れます．これらは属性情報付きの地図上の図形情報を表すためのデータ形式，"SHAPE形式"だそうです．この中に駅別乗降客数の情報も入っているということで，これを取り出していきます（GeoJSONもありましたが今回は使わなかった）．

:::note info
**SHAPE形式をPythonで読み込む**
SHAPE形式をPythonで読み込むには`shapefile`パッケージが便利です．
:::


```python
import shapefile
sf = shapefile.Reader("./S12-23_GML/utf8/S12-23_NumberOfPassengers.shp")
```
`sf.shape`で図形情報，`sf.record`で属性情報を取得できます．

図形情報としては，駅の端と端を結んだ線分が収納されているみたいです：
```python
sf.shape(0)
```
```text
Shape #0: POLYLINE
```
```python
sf.shape(0).points
```
```text
[(130.63035, 31.25405), (130.62985, 31.25459)]
```
属性情報としては，駅情報の諸々と駅別乗降客数が記載されているようです：
```python
sf.record(0)
```
```text
Record #0: ['二月田', '010112', '010112', '九州旅客鉄道', '指宿枕崎線', 11, 2, 1, 3, '', 0, 1, 3, '', 0, 1, 3, '', 0, 1, 3, '', 0, 1, 3, '', 0, 1, 3, '', 0, 1, 3, '', 0, 1, 1, '', 658, 1, 1, '', 690, 1, 1, '', 318, 1, 1, '', 305, 1, 1, '', 622]
```
で，この各カラムは`sf.fields`に記載されています：
<details><summary>詳細</summary>

```text
[('DeletionFlag', 'C', 1, 0),
 ['S12_001', 'C', 254, 0],
 ['S12_001c', 'C', 254, 0],
 ['S12_001g', 'C', 254, 0],
 ['S12_002', 'C', 254, 0],
 ['S12_003', 'C', 254, 0],
 ['S12_004', 'N', 11, 0],
 ['S12_005', 'N', 11, 0],
 ['S12_006', 'N', 11, 0],
 ['S12_007', 'N', 11, 0],
 ['S12_008', 'C', 254, 0],
 ['S12_009', 'N', 11, 0],
 ['S12_010', 'N', 11, 0],
 ['S12_011', 'N', 11, 0],
 ['S12_012', 'C', 254, 0],
 ['S12_013', 'N', 11, 0],
 ['S12_014', 'N', 11, 0],
 ['S12_015', 'N', 11, 0],
 ['S12_016', 'C', 254, 0],
 ['S12_017', 'N', 11, 0],
 ['S12_018', 'N', 11, 0],
 ['S12_019', 'N', 11, 0],
 ['S12_020', 'C', 254, 0],
 ['S12_021', 'N', 11, 0],
 ['S12_022', 'N', 11, 0],
 ['S12_023', 'N', 11, 0],
 ['S12_024', 'C', 254, 0],
 ['S12_025', 'N', 11, 0],
 ['S12_026', 'N', 11, 0],
 ['S12_027', 'N', 11, 0],
 ['S12_028', 'C', 254, 0],
 ['S12_029', 'N', 11, 0],
 ['S12_030', 'N', 11, 0],
 ['S12_031', 'N', 11, 0],
 ['S12_032', 'C', 254, 0],
 ['S12_033', 'N', 11, 0],
 ['S12_034', 'N', 11, 0],
 ['S12_035', 'N', 11, 0],
 ['S12_036', 'C', 254, 0],
 ['S12_037', 'N', 11, 0],
 ['S12_038', 'N', 11, 0],
 ['S12_039', 'N', 11, 0],
 ['S12_040', 'C', 254, 0],
 ['S12_041', 'N', 11, 0],
 ['S12_042', 'N', 11, 0],
 ['S12_043', 'N', 11, 0],
 ['S12_044', 'C', 254, 0],
 ['S12_045', 'N', 11, 0],
 ['S12_046', 'N', 11, 0],
 ['S12_047', 'N', 11, 0],
 ['S12_048', 'C', 254, 0],
 ['S12_049', 'N', 11, 0],
 ['S12_050', 'N', 11, 0],
 ['S12_051', 'N', 11, 0],
 ['S12_052', 'C', 254, 0],
 ['S12_053', 'N', 11, 0]]
```
</details>

0要素目のDeletionFlagとかいうのは無視して，1要素目以降にカラムの情報が記載されています．各要素の0要素目がカラム名だそうですが，よく分からんコード形式になっています：
```python
fieldname_list = [field[0] for field in sf.fields[1:]]
fieldname_list
```
```text
['S12_001',
 'S12_001c',
 'S12_001g',
 'S12_002',
 'S12_003',
 'S12_004',
 'S12_005',
 ...
 ```

 カラムの並びについては，国交省の[ダウンロードページ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-S12-2022.html)の中に書かれています：
 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3116310/f78e45d3-f99c-4764-b6ae-0987e2761ece.png)

分かりやすさのため，カラム名の対応表を作っておきます．今回は駅名，企業名，路線名，各年の乗降客数，状況客数データの有無を対応づけておきます．
```python
FIELDNAME_TO_COLNAME = {
    "S12_001": "StationName",
    "S12_002": "CompanyName",
    "S12_003": "LineName",
    "S12_039": "NumPassengersDataStatus19",
    "S12_041": "NumPassengers19",
    "S12_043": "NumPassengersDataStatus20",
    "S12_045": "NumPassengers20",
    "S12_047": "NumPassengersDataStatus21",
    "S12_049": "NumPassengers21",
    "S12_051": "NumPassengersDataStatus22",
    "S12_053": "NumPassengers22",
}
FIELDIDX_TO_COLNAME = {
    fieldname_list.index(fn): cn
    for fn, cn in FIELDNAME_TO_COLNAME.items()
}
```

続いて，これらに挙げたカラム名と駅の緯度経度をPandasのDataFrameに転記したいと思います：

:::note info
`DataFrame`のコンストラクタにはGeneratorを渡すことができます．
:::

```python
def create_row(record, shape) -> dict:
    row = {}
    for i, cn in FIELDIDX_TO_COLNAME.items():
        row[cn] = record[i]
    row["Lon1"] = shape.points[0][0]
    row["Lat1"] = shape.points[0][1]
    row["Lon2"] = shape.points[1][0]
    row["Lat2"] = shape.points[1][1]
    return row

rows_gen = (
    create_row(record, shape)
    for record, shape
    in tqdm(
        zip(sf.iterRecords(), sf.iterShapes())
    )
)
df = pd.DataFrame(rows_gen)
```

CSVに保存しようと思います．

:::note info
**CSV出力の際のカラム順序**
`DataFrame.to_csv`メソッドでは，`columns`引数を加えることによってカラムの順序を指定できます．
:::

```python
column_order = """StationName
CompanyName
LineName
Lat1
Lon1
Lat2
Lon2
NumPassengersDataStatus19
NumPassengers19
NumPassengersDataStatus20
NumPassengers20
NumPassengersDataStatus21
NumPassengers21
NumPassengersDataStatus22
NumPassengers22""".splitlines()
df.to_csv("./StationPassengerData.csv", columns=column_order)

```

CSVファイルの内容はこんな感じになります：
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3116310/30093cbd-e635-40d3-8965-c2ce96c4c277.png)

駅は全部で10,500もありました！！

## 2. デザインの準備
### 2.1. シンプルな指数表記
GoogleEarthでは，各駅の位置の上に乗降人数を指数表記で表示しようと思います．なにせ1万個も駅があるので，1駅あたりの表示はなるべくシンプルにしたいです．

Pythonのformatにも指数表記はありますが，例えば`3e+00`のようにちょっと冗長な表記になってしまうんです．次のような極力シンプルな指数表記を自作しようと思います：

元の数 | 目指す指数表記
--- | ---
3 | 3e0
24 | 2e1
98 | 1e2
103 | 1e2
9382 | 9e3

最も上の位を丸め，どんな数もよりシンプルに表すようにします．

で，これが実装です：
```python
def divide_base_and_exp(x: int) -> Tuple[int, int]:
    if x == 0:
        return 0, 0
    exp = int(math.floor(math.log10(x)))
    base = round(x, -exp)//(10**exp)
    if base == 10:
        exp += 1
        base = 1
    return base, exp


def simple_sci(x: int) -> str:
    base, exp = divide_base_and_exp(x)
    return f"{base}e{exp}"
```

### 2.2. カラースケール
GoogleEarthでは，乗降客数に応じて色を変えるようにしたいです．そのために，Matplotlibのカラーマップを使おうと思います．
例えば，次の例では，2022年における乗降客数をSpringカラーマップで表示することを考え，色を用意しています．
```py
cm = colormaps["spring"]
num_passengers = df.NumPassengers22.to_numpy()
point_colors = cm(num_passengers / num_passengers.max()*cm.N)
point_colors_int = (point_colors*255).astype(int)
kml_colors = [simplekml.Color.rgb(r, g, b) for r, g, b, a in point_colors_int]
```

:::note info
KML形式では色データの順番はRGBではありません．`simplekml.Color.rgb`を使用して，RGBからKML形式の色データを作ります．
:::


## 3. KMLを作る

```python
kml = simplekml.Kml()
for (i, record), num_passenger, color in tqdm(zip(df.iterrows(), num_passengers, kml_colors)):
    if record.NumPassengersDataStatus22 != 1:
        name = "."
        if num_passenger == 0:
            continue
    else:
        name = simple_sci(num_passenger)

    point = kml.newpoint(name=name)
    point.coords = [(
        (record.Lon1 + record.Lon2)/2,
        (record.Lat1 + record.Lat2)/2
    )]
    point.style.labelstyle.color = color
    point.style.labelstyle.scale = 1
    point.style.iconstyle.icon.href = ""

kml.save("StationPassengersMap.kml")
```

## 4. 完成
出てきたKMLファイルをGoogleEarthに読み込ませると，いい感じに表示してくれます．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3116310/bae12e03-9aa8-412c-a65c-d4eef348fd0d.png)