pythonで地域メッシュごとのポリゴンデータからGeoTIFFを作成する：国土数値情報平年メッシュデータを例に - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

TL;DR

地域メッシュポリゴンデータからGeoTIFFを作る手順を解説
pythonのrioxarrayを使うと簡単
テストケースで作成した平均年間雨量のデータが興味深い。

結果だけ参照されたい方は「国土数値情報：平均年間雨量メッシュデータ」まで読み飛ばして下さい。

モチベーション
地域メッシュとは？
元shapefileファイルのデータ構造
GeoTIFF用二次元配列の作成
メッシュセンターの緯度経度の設定
GeoTIFFの作成
GeoTIFFの結合
- VRTの作成
- 1つのGeoTIFFを作成
国土数値情報：平均年間雨量メッシュデータ
GitHub

モチベーション

国土数値情報 | 平年値メッシュデータが、1次地域メッシュごとのshapefileでポリゴンデータというかなり扱いにくい形式だったので、この手のデータの近年の主流であるGeoTIFFに変換してみました。

地域メッシュとは？

過去記事欠損したグリッド標高テキストデータからgeotiffを作成するpythonスクリプト - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning で解説した航測測量データの図郭は平面直角座標系（XY）で定義されますが、地域メッシュは地理座標系（緯度経度）で下表のとおりに定義されます。（出典はこちら）左下を起点に番号が振られています。

参照：統計局ホームページ/地域メッシュ統計の地域メッシュ統計の特質・沿革

元shapefileファイルのデータ構造

shapefileファイルのデータ構造の主な特徴は次のとおりです。（詳細はこちらを参照）

1次メッシュごとに153個のファイル作成されており、各ファイルに3次メッシュごとのポリゴンが格納されています。
そのため、各ファイルのポリゴン数は80X80になりますが、水域はデータが作成されていないため、ファイルごとにポリゴン数が異なります。
各ポリゴンの属性に3次メッシュコードが格納されています。

GeoTIFF用二次元配列の作成

80x80の二次元配列を作成して3次メッシュコードから値を入力します。型はuint16（範囲は0～65535）を使用しますので、nodataを65000にしています。

こちらのG02-12_5337-jgd_GMLのファイルから平均年間雨量を雨量のデータからGeoTIFFを作ります。

import geopandas as gpd
import numpy as np

fp = "G02-12_5337-jgd.shp"
gdf = gpd.read_file(fp)

ind = gdf["G02_001"].values # 3rd mesh code
val = gdf["G02_014"].values # Total annual rainfall

# size 1st order Chiki-Mesh
x = np.full((80,80), int(65000), dtype='uint16') # y, x

for dp, v in zip(ind, val):
    # yindex
    iy = 10 * int(dp[4]) + int(dp[6])
    # yindex
    ix = 10 * int(dp[5]) + int(dp[7])
    x[iy,ix] = int(v)

メッシュセンターの緯度経度の設定

前述のとおり、メッシュコードからポリゴン左下の座標を計算できます。値はメッシュのセンターで定義するため緯度経度方向ともに半メッシュずらします。

dp = ind[0] #3次メッシュコード
lonp = float(100) + float(dp[2:4])  + float(0.5*45/3600)
latp = float(dp[:2]) * float(40/60) + float(0.5*30/3600)

lon = lonp + np.arange(80)*45/3600
lat = latp + np.arange(80)*30/3600

GeoTIFFの作成

GeoTIFFの作成はGDALを使いますが、（現時点では）Pythonicな記述ができないため、GDALがwrapされているrioxarray、fiona等を使ったほうがスマートです。今回はrioxarrayを使います。

GeoTIFFの圧縮は、画像処理のど素人がまとめたpythonによるTIFF圧縮形式の比較 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning を参照頂きたいです。

import xarray as xr
import rioxarray

# y is North to South
ds = xr.Dataset({
                 'value': (['y','x'], x[::-1])                 
                }
                , coords={
                            'x': lon
                          , 'y': lat[::-1]
                         }
                    , attrs={'crs':'+init=epsg:' + str(4326)}
                )

dsp = ds['value']
d = dsp.rio.write_crs('epsg:4326', inplace=True)
out = dsp.rio.to_raster(f'geotiff\\{dp[:4]}.tif') #, compress='zstd')

GeoTIFFの結合

上記の手順で153個全てのshapefileをGeoTIFFに変換した後、取り扱いやすいように一つのファイルにまとめます。

方法は以下の2種類が考えられます。

VRTの作成

お馴染みのVRTを使う方法です。個々のファイルが重い場合はこの方法が良いです。

from osgeo import gdal 

# gdal 2.1以降
my_vrt = gdal.BuildVRT('geotiff/output.vrt', glob.glob( 'geotiff/*.tif'), VRTNodata=65000, srcNodata=65000)
my_vrt = None

1つのGeoTIFFを作成

今回は個々のファイルが軽いので一つのファイルにします。上記で作成したVRTから変換しています。

ds = xr.open_rasterio('geotiff/output.vrt')

d = ds.rio.write_crs('epsg:4326', inplace=True)
out = ds.rio.to_raster('geotiff\\AverageAnnualRainfall2010.tif', compress='zstd')

試してはいませんが、gdal_mergeでも出来るようです。（参照：python - Stack geotiff images while retaining individual bands - Stack Overflow）