gdal_gridでgeojsonをバイナリに変換する


gdalフレームワークの`gdal_grid`でgeojsonファイルをバイナリデータに変換する方法を紹介します。

## gdalのインストール
まずは`gdal`をインストールしておきます。
`macbrew`が入っているmacであれば以下の通り。

```
brew install gdal
```

ubuntuであれば以下の通り。

```
sudo add-apt-repository ppa:ubuntugis/ppa
sudo apt-get update
sudo apt-get install gdal-bin
```

インストールできているかの確認。
`gdal_grid -h`と入力したら、ヘルプが表示されればインストール成功です。

```
gdal_grid -h
>>> Usage: gdal_grid [--help-general] [--formats]
    [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/ ...
```

## gdal_gridとは
`gdal_grid`は散布状に存在している点データを格子状のデータに変換してくれるコマンドです。
元のデータに抜け落ちがあったり、そもそも格子状に散布されていない場合でも、適切な補完を選択することで格子状に再グリッドしてくれます。
![スクリーンショット 2018-07-20 17.57.35.png](https://qiita-image-store.s3.amazonaws.com/0/61002/139df38d-cfe1-10c7-befb-f4815aebb6d1.png "スクリーンショット 2018-07-20 17.57.35.png")


オプションは以下の通り。

>
gdal_grid [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
          CInt16/CInt32/CFloat32/CFloat64}]
          [-of format] [-co "NAME=VALUE"]
          [-zfield field_name] [-z_increase increase_value] [-z_multiply multiply_value]
          [-a_srs srs_def] [-spat xmin ymin xmax ymax]
          [-clipsrc <xmin ymin xmax ymax>|WKT|datasource|spat_extent]
          [-clipsrcsql sql_statement] [-clipsrclayer layer]
          [-clipsrcwhere expression]
          [-l layername]* [-where expression] [-sql select_statement]
          [-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]
          [-a algorithm[:parameter1=value1]*] [-q]
          <src_datasource> <dst_filename>
公式ドキュメント：[http://www.gdal.org/gdal_grid.html](http://www.gdal.org/gdal_grid.html)

## 変換方法
ここでは例として、国土地理院の5m解像度の標高geojsonデータ（ダウンロード方法などの説明：[https://github.com/gsi-cyberjapan/experimental_dem](https://github.com/gsi-cyberjapan/experimental_dem)）をバイナリに変換します。

以下のコマンドでgeojsonファイルを変換できます。

```
gdal_grid -ot Float32 -of ENVI -l OGRGeoJSON -zfield alti -a nearest -txe $lonmin $lonmax -tye $latmin $latmax -outsize 600 800 geojson_name bin_name
```

* `-ot Float32`: 出力されるバイナリデータ内の１データごとの保存形式。float32形式を指定。
* `-of ENVI`：出力されるファイル形式を指定。ピュアなバイナリデータのファイルが欲しいのでここではENVIを指定。
* `-l OGRGeoJSON`：元データのどのレイヤを参照するかを指定。`OGRGeoJSON`もしくはファイル名を指定する。詳しくは次節で後述。
* `-zfield alti`：元データのどの値を出力するかを指定。本タスクの場合、元データは以下のようになっているので、`alti`を指定。

```
{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "geometry": {
                "type": "Point",
                "coordinates": [
                    137.6285,
                    36.5613888888889
                ]
            },
            "properties": {
                "class": "DEMPt",
                "type": "\u305d\u306e\u4ed6",
                "alti": 2250.9
            }
        }, ...
```

* `-a nearest`：補完アルゴリズム。いろいろ選べるが、今回は最近傍法で補完するのでnearestを指定。
* `-txe $lonmin $lonmax -tye $latmin $latmax`：出力する範囲の経度、緯度を指定する。ここでは`$lonmin`とかにすでに値が入っている想定。
* `-outsize 600 800`：出力バイナリデータのピクセル数。横ピクセル数、縦ピクセル数の順に指定。
* `geojson_name`：入力geojsonデータ名。
* `bin_name`：出力バイナリファイル名。


## `-l`レイヤ名の指定
以下、`gdal_grid`で指定した`-l OGRGeoJSON`オプションに関してです。
geoJSONにはレイヤはなかったのに、何を指定すればいいんだ？？？となるかもしれません。

gdalでは実はgeojsonを読み込むときに、１ファイル１レイヤのOGRLayer objectとして取り込んでいます。
そのため、そのレイヤの名前をしてあげる必要が出てくるわけです。（１つのレイヤしかないのならばどうにかしてくれればとも思うのですが...）

そのレイヤの命名規則は以下の通りです。

```
1. "FeatureCollection"レベルに"name"という属性があれば、それがレイヤ名。
2. 1を満たさなければ、ファイル名がレイヤ名になる。ただし、ファイル名が特殊（"."があったりするとダメらしい）だとダメ。
3. 1も2も満たさなければ、レイヤ名は"OGRGeoJSON"になる。
```
上記の例ではいろいろとあって、レイヤ名は"OGRGeoJSON"になっているようでした。

参考：[http://www.gdal.org/drv_geojson.html](http://www.gdal.org/drv_geojson.html)


## 国土地理院　標高データ変換コード
`gdal_grid`を使って、国土地理院の標高データ（DEM）をバイナリに変換するコードを作成しました。
欲しい範囲の座標さえ入力すれば、必要データのダウンロードから整形、バイナリデータに変換、確認用画像生成まで一気にやってくれます。
コードをgithubにおいておいたので、参考にしたり使ったりしてください。
[https://github.com/mski-iksm/kokudoDEM2bin](https://github.com/mski-iksm/kokudoDEM2bin)

```
git clone git@github.com:mski-iksm/kokudoDEM2bin.git
cd kokudoDEM2bin
python kokudoDEM2bin.py -l 35.519373 139.741287 35.581384 139.854240
```