見出し画像

QGISを用いた土地利用被覆度の計算

みつばちさん

本記事では、地図データを扱うことができるソフトQGISを用いて、JAXAが提供する高解像度土地利用土地被覆図のデータから、調査地点の半径3 km圏内の土地利用状況を算出する方法を解説します。
最終的には以下のような図をQGISで表示して、土地利用被覆度を計算するのがゴールです。

土地利用被覆度（Ratio）を求めた結果が以下の表です。
東京タワーから半径3 kmの範囲は、建物（Built-Up）が81.63%を占めており、落葉針葉樹林（DNF）がほとんどない（0.02%）などがわかります。


本記事の実行環境

OS：macOS Sonoma version 14.6.1

事前準備

  1. QGISのインストール

  2. JAXAの高解像度土地利用土地被覆図のダウンロード

  3. 調査地点の座標を得る（Google map）

1. QGISのインストール

この方法では、QGISという無料ソフトを利用して地図データを処理します。
以下のサイトから自分の利用環境（Windows, Mac, Linux…）にあったバージョンをダウンロードします。

2. JAXAの地図データのダウンロード

JAXAの高解像度土地利用土地被覆図を利用するために、以下のサイトで利用登録します。

特に理由がなければ、最新のバージョンをダウンロードすればいいです。
ダウンロードしたZipファイルを解凍すると、緯度(N)と経度(E)で分割されたGeoTIFF形式のデータが大量に確認できると思います。自分が利用したいデータのみ今後は利用します。（例：東京なら"N35E139"のデータ）

3. Google mapから調査地の座標を取得

今回は東京タワーを調査地点として、その半径3 kmの土地利用状況を調べたいと思います。
Google mapを開き、東京タワーを検索し、右クリックすると座標(緯度経度)が表示されます。この座標をコピーして、CSV ファイルにまとめます。

IDをtokyo tower、先ほどの座標のNに緯度、Eに経度をコピペします。

CSV形式で座標データを保存します。
これで事前準備は完了！

QGISに地図データを追加する

地図データを追加する

QGISを立ち上げたら、左下にあるレイヤに先ほどJAXAからダウンロードしたGeoTIFF形式の地図データをドロップします。

プロジェクトのCRSの設定を変更する

次にプロジェクトの座標参照(CRS)をメートル法に変更します。右下にある「EPSG:xxxx」をクリックします。すると、プロジェクトのプロパティ座標参照系(CRS)が開きます。ここで座標参照系（CRS）を「JGD2011 / Japan Plane Rectangular CS Ⅸ, EPSG : 6677」に設定します。

これでこのプロジェクトの設定が変更できました。（一時的に地図データが見られなくなるかもしれません）
この設定により、この後で必要になるメートル法による計算ができるようになります。

地図データのCRS設定を変更する

続いて、先ほどレイヤに追加した地図データのCRSを変更します。
レイヤの地図データで右クリックし「エクスポート」から「名前をつけて保存」を選択します。ファイル名を設定しますが、パスまで入力が必要なので、記入欄の右の「•••」を押して、保存したいフォルダを選択してファイル名を入力した方が楽です。で、座標参照系（CRS）を先ほどと同じ設定「EPSG : 6677 JGD2011」に設定し保存します。

保存できるとレイヤパネルに、CRSが変更された新しく地図データが追加されます。
これで、プロジェクトと地図データのCRSが同じ設定になりました。

QGISに調査地点のデータを追加する

CSVファイルをレイヤに追加する

「レイヤ」タブ内にある「レイヤを追加」から「CSV テキストレイヤを追加」を選択する。

一番上にあるファイル名の右端にある「•••」をクリックして、先ほど作成した調査地のCSV ファイルを選択する。
ジオメトリ定義にあるジオメトリのCRSを「EPSG : 4326 WGS 84」に設定します。東京タワーの緯度と経度の座標（地理座標系）で取得したので、それに合わせたCRSを設定し、東京タワーの座標データを追加します。

追加すると地図の上に調査地点が⭕️で表示される。

調査地点（tokyo_tower）のレイヤをshape形式で保存する

レイヤパネルに追加された「tokyo_tower」で右クリックし、「エクスポート」から「新規ファイルに地物を保存」に進みます。

形式を「ESRI Shapefile」に設定し、ここでCRSを「EPSG:6677(プロジェクトと同じ)」に変更して保存します。
そのまま「OK」を押すとレイヤに新しい調査地のレイヤが追加されます。

レイヤパネルに新しい「tokyo_tower_location」のレイヤが追加されたら、先に追加した方は削除します。

地図データをベクタ化する

地図データをベクタ化します。「ラスタ」タブ内にある「変換」から「ラスタをベクタ化(polygonize)」に進みます。

入力レイヤに地図データを選択し、実行をします。
ラスタ化には10分ぐらいの時間がかかります。

ラスタ化が完了するとレイヤパネルに「ベクタ化」という新しい地図データが出力されます。
正しく表示されていたら、元の地図データは削除します。
新しく作成された「ベクタ化」レイヤは一時ファイルなので、保存する前にQGISを閉じるとデータが消えてしまいます。なので「ベクタ化」レイヤで右クリックして、Shape形式で保存します。

バッファを作成する

「バッファ」機能を用いることで、調査地点から半径3 kmの範囲のレイヤを作成することができます。
「ベクタ」タブ内の「空間演算ツール」から「バッファ(buffer)」に進みます。

入力レイヤに調査地点のshapeファイルを設定し、距離を3000メートルに設定し、実行します。（この処理もしばらく時間がかかります）

バッファが完了したらレイヤに新しい「出力レイヤ」という名前のレイヤが追加されます。このレイヤも一時レイヤなので、保存しておきます。

指定した距離（3000 m）に応じたレイヤが地図の上に表示されたら成功です。

交差(interaction)の実行

「ベクタ」タブ内にある「空間演算ツール」から「交差(intersect)」に進む。
入力レイヤには、「ベクタ化」を設定し、オーバーレイヤには、先ほどバッファを実行した際に追加された出力レイヤを設定する。

実行すると、新しく「交差 (intersect)」レイヤが追加されます。このレイヤで右クリックし「属性テーブルを開く」を押します。

属性テーブルが開いたら「フィールド計算機を開く」を押します。
ここで、新しいフィールド（列）を追加することができます。
フィールドの名前を「Area」にし、フィールド型を「小数点付き数値」に設定します。次に、ジオメトリ内にある「$area」関数をダブルクリックして式として追加し実行します。

これで、東京タワーから半径3 km圏内の各土地利用方法（各DN）の面積を計算した列（フィールド）が新しく属性テーブルに追加されます。

交差レイヤの保存

続いて、レイヤパネルの交差レイヤで右クリックし、「エクスポート」から「新規ファイルに地物を保存」に進みます。

形式を「カンマで区切られた値 CSV」に設定し、ファイル名を指定して保存します。

googleの衛星写真を追加

最後はgoogleの衛生写真を追加して綺麗に仕上げます。
ブラウザパネルにある「XYZ Tiles」で右クリックし、「新規接続」に進みます。

名前「google satellite map」（好きな名前でOK）を入力し、URL欄に以下のURLを入力しOKします。

https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}

XYZ Tilesに新しくgoogle satellite mapが追加されます。これをダブルクリックまたはレイヤパネルにドラッグすれば地図パネルに衛生写真が追加されます。出力レイヤの円を薄くしたい場合は、右クリックからプロパティへ進み、不透明度を下げる（下の図は70%）ことで、地図も透けて見えるようになります。

この地図を保存したい場合は、メニューの「プロジェクト」タブ内にある「インポートとエクスポート」に進み保存することができます。

これでQGISでの処理は終わりです👏

Pythonで土地利用被覆度を求める

交差の実行により、東京タワーから半径3 kmの各土地利用方法の面積が得られました。次はこのデータから、土地利用被覆度（比率）を求めたいと思います。エクセルなどの表計算ソフトで十分可能な計算ですが、今回はPython（jupyter Notebook）で計算してみます。

jupyter Notebook

jupyter NotebookのHome画面の右上にある「Upload」から先ほど保存した属性テーブルのデータをアップロードします。
「New」から新しいNotebookを作成し、以下のコードを入力します。
実行すれば結果が表示されるはずです。

import pandas as pd

# CSVファイルを読み込む
file_name = "table.csv" # 保存した属性テーブルのファイル
df = pd.read_csv(file_name)

# Area列の合計を計算
total_area = df['Area'].sum()

# DNの値を変換（DNに対応する土地利用方法は地図データをダウンロードした際のJAXAのフォルダ"README"にある）
replacements = {
    1: 'Water bodies',
    2: 'Built-up',
    3: 'Paddy field',
    4: 'Cropland',
    5: 'Grassland',
    6: 'DBF (deciduous broad-leaf forest)',
    7: 'DNF (deciduous needle-leaf forest)',
    8: 'EBF (evergreen broad-leaf forest)',
    9: 'ENF (evergreen needle-leaf forest)',
    10: 'Bare',
    11: 'Bamboo forest',
    12: 'Solar panel',
    13: 'Wetland',
    14: 'Greenhouse'
}
df = df.replace(replacements)

# 使わない列を削除
df = df.drop(['fid', 'N', 'E'], axis=1)

# 土地利用方法（DN）の合計面積を計算
DN_totals = df.groupby('DN')['Area'].sum()

# Area合計面積を計算
total_area = df['Area'].sum()

# 比率を計算
ratios = (DN_totals / total_area) * 100

# インデックスを列としてデータフレームに変換
ratios = ratios.reset_index()

# 列の名称を追加
ratios.columns = ['Category', 'Ratio']

# 比率を小数点以下2桁にフォーマット
ratios['Ratio'] = ratios['Ratio'].map(lambda x: round(x, 2))

print(ratios)

各DNの値が何を示しているかは、JAXAからダウンロードしたフォルダ内にある「 README」ファイルに示されています。
データ"ratios"を表にすると以下のような表ができあがります。

よくあったエラー

レイヤが表示されない

QGISではCRSの設定が合っていないとレイヤのデータが地図パネルに正しく反映されません。レイヤが表示されていない場合は、プロジェクトのCRSとレイヤのCRSが一致しているかを確認します。
あとは、「ベクタ化」や「バッファ」などの処理で出力されたレイヤは「保存」や「エクスポート」をしてから、再び新しいレイヤとして追加する。

フリーズする

QGISの各処理は重たい作業になるので、PCにかなり負担がかかります。したがって、たまにフリーズしてしまいます。そんな時は、強制終了するしかないので、QGIS閉じると消去されてしまう一時レイヤはすぐに保存しておく。

proj.dbが見つからない

proj.dbはCRSに関するファイルのようです、このファイルがないと上手くデータを変換することができません。「proj.dbが見つからない」というエラーが生じた場合

  1. PC内にproj.dbファイルが存在するか確認する

  2. proj.dbファイルがあれば、パスを設定し直す

  3. proj.dbファイルがなければ、ダウンロードする


いいなと思ったら応援しよう！