QGISで緯度経度を指定したグリッド（メッシュ）作成と使用例

QGIS　バージョン3.40

1.はじめに

今回は、QGISで緯度経度を指定したグリッドを作成していきたい思います。使用用途としては、メッシュデータの作成を想定しています。国土数値情報等でダウンロードできるデータは、比較的に使い易い形になっているのですが、中にはうまく表示できないものもあるので、力技で手間がかかるのですが、ご参考になれば幸いです。今回は例として、最近標高改定前後のジオイド2011とジオイド2024をQGIS上で表示してみたいと思います。

QGISには、デフォルトの機能としてグリッドを作成する機能があるのですが、添付した画像のように、緯度経度ではなくXY座標を指定します。また、各方向の間隔も緯度経度から距離換算した場合、切りのいい数字にならないため、デフォルトの機能では緯度経度を指定したグリッドの作成が難しいです。

2.グリッド作成

QGISを起動しプラグイン→Pythonコンソールから下記のコードを実行します。

# 北西を起点に東方向
import qgis.utils
from qgis.core import (
    QgsRectangle,
    QgsVectorLayer,
    QgsVectorDataProvider,
    QgsField,
    QgsFields,
    QgsGeometry,
    QgsFeature,
    QgsPointXY,
    QgsProject,
)
from PyQt5.QtCore import QVariant

# グリッドを作成する範囲を指定 (座標系はEPSG:4326)
xmin, xmax, ymin, ymax = 119.9875, 160.0125, 14.9916665, 50.0083335

# X方向、Y方向の分割数
nx, ny = 1601, 2101
width = (xmax - xmin) / nx
height = (ymax - ymin) / ny

layer_name = "grid"
crs = QgsProject.instance().crs()
grid = QgsVectorLayer("Polygon?crs=EPSG:4326", layer_name, "memory")
provider = grid.dataProvider()

fields = QgsFields()
fields.append(QgsField("ID", QVariant.Int))
provider.addAttributes(fields)
grid.updateFields()

feature_list = []
feature_id = 1

for j in range(ny-1, -1, -1):
    for i in range(nx):
        xmin_, ymin_ = xmin + i * width, ymin + j * height
        xmax_, ymax_ = xmin + (i + 1) * width, ymin + (j + 1) * height
        
        coords = [
            QgsPointXY(xmin_, ymin_),
            QgsPointXY(xmax_, ymin_),
            QgsPointXY(xmax_, ymax_),
            QgsPointXY(xmin_, ymax_),
            QgsPointXY(xmin_, ymin_)
        ]
        
        polygon = QgsGeometry.fromPolygonXY([coords])
        feature = QgsFeature()
        feature.setGeometry(polygon)
        feature.setAttributes([feature_id])
        feature_list.append(feature)
        feature_id += 1

batch_size = 10000
for i in range(0, len(feature_list), batch_size):
    provider.addFeatures(feature_list[i:i+batch_size])

QgsProject.instance().addMapLayer(grid)
print(f"グリッド作成完了: {feature_id-1}個のセルを作成しました")

作成したグリッド（メッシュ）にジオイド2024のデータを結合して表示します。

以上がグリッドの作成になります。以降は、作成したグリッドの使用例になります。

3.データ準備

国土地理院ジオイドモデル入手より、ジオイド2011、2024をダウンロードします。ジオイド2024は基盤地図情報からダウンロードするので、国土地理院のアカウント登録が必要になります。

ジオイド2011

ジオイド2024

上の画像のようにジオイド2011とジオイド2024でデータの形や拡張子が異なります。また、ジオイド2024は北西を起点に東方向にデータが並んでいますが、ジオイド2011は、南西を起点に北方向にデータが並んでいます。さらに、ジオイド2011は、１行が43行に分けられているため、データの加工が必須になります。少しややこしいですが、上記の情報がわかっていればそこまで難しくないので、データ加工は割愛します。
それぞれ１列に整列し、IDを振ってCSVにします。

データの準備はこれにて完了です。IDを振ったのは、グリッドのIDとリンクさせて数値を入れるためです。ちなみに、ジオイド2024は用途にもよりますがほぼ加工せずに使用できます。

グリッド作成のために必要なパラーメータは下記になります。

  • 範囲（緯度経度）
  • 分割数

今回の場合、下記のようにとなります。

ジオイド2011
緯度　19.9916665°～　50.0083335°分割数　1801
経度　119.9875°～　150.0125°分割数　1201

ジオイド2024
緯度　14.9916665°～　50.0083335°　分割数　2101
経度　119.9875°～　160.0125°分割数　1601

【上記の緯度経度と分割数の算出方法】
ジオイド2011を例に分割数の考え方と緯度経度の算出方法を解説します。
gsigeo2011_ver2_2.ascの１行目は解説PDFに記載がある通り下記のようになっています。

緯度経度の間隔が馴染みない数字になっていますが、これは度分秒表示を度の十進法に直した数字になります。

緯度の間隔　0°01'00"＝ 1/60°≒ 0.016667
経度の間隔　0°01'30"＝ 1/40°＝ 0.025000

分割数と緯線経線の個数についてですが、普通に考えると分割数は緯線経線の本数に比べ1少なくなるのですが、データ数を確認すると分割数と緯線経線の本数が同じになっています。このことから、グリッドの中心を軸とした考え方をしており、メッシュを配置した際、メッシュの端は緯度経度の間隔の半分外側にはみ出ることになると思われます。一応調べましたが、確証できるものがなかったので、不安な方は国土地理院に問い合わせをしてみてください。
分割数と緯度経度の間隔から緯度経度の終点を求めると下記のようになります。

緯度
緯度の間隔　1/60°×　分割数　1801　＝　30.016667°
緯度起点　20°-　（0.016667°/2）＝　19.9916665°
緯度終点　19.9916665°+　30.016667°＝　50.0083335°
　　　　
経度
経度の間隔　1/40°×　分割数　1201　＝　30.025°
経度起点　120°-　（0.025°/2）＝　119.9875°
経度終点　119.9875°+　30.025°＝　150.0125°

4.グリッド作成

前述したパラーメータを元にグリッド作成します。コード内の下記の箇所を編集することで、グリッドをどこからどの方向に作成するか選択できます。

#該当箇所
for j in range(ny-1, -1, -1):
    for i in range(nx):
西から東
for i in range(nx):
東から西
for i in range(nx-1, -1, -1):
南から北
for j in range(ny):
北から南
for j in range(ny-1, -1, -1):

南西を起点に北方向
for i in range(nx):
    for j in range(ny):
北東を起点に西方向
for j in range(ny-1, -1, -1):
    for i in range(nx-1, -1, -1):

ジオイド2011を例にグリッドを作成します。

5.CSVデータとの結合

作成したグリッドと準備したCSVデータを結合します。レイヤ→レイヤを追加→CSVテキストレイヤを追加からCSVデータを読み込みます。この時、ジオイド定義をジオメトリなし（属性テーブルのみのテーブル）にします。

gridレイヤのプロパティを開き、テーブル結合を選択します。「+」を押してテーブル結合をします。軽くテーブル結合について触れておくのですが、機能としてはExcelのVlookupに近いです。共通のフィールドを参照し、対応する他のフィールドを結合できます。CSVデータにIDを振ったのは、メッシュのIDと共通のフィールドを作成する必要があったからです。

結合するレイヤとフィールドを選択します。

gridの属性テーブルを開くと下のようになります。

シンボリジを編集していい感じにします。

6.おわりに

今回は、緯度経度を指定したグリッド作成してみました。グリッドの間隔が割り切れない数値だと正確に作成できないのが難点ですが、多くのメッシュを基準にしたデータに対しての解決方法になると思います。

最後にジオイド2011とジオイド2024を比較したものを載せておきます（ちょっと失敗しています）。
ジオイド2011とジオイド2024は、緯度経度の範囲やデータの並び順が違いますが、メッシュに割り当てることで、視覚的わかりやすく、同一メッシュで差を比較することもできます。また、フィールド計算機を使うことで同一地点のジオイドの変動量を算出できます。メッシュの結合は、和集合でしているのですが、かなりスペックの高いPCでも1時間ほどかかるので注意が必要です。

下は、ジオイドの変動量を示したものになります。

コメント