J-EGG500の水深データをGeoTiffに変換しました(補間がいまいち)

秋田港周辺の海底地形を作成したかったので、まず"海底地形データ フリー"で検索したところ

GMTTipsPage

のリンクがありJ-EGG500というデータが日本周辺では定番のようです。早速、日本海洋データセンターのサイトに行ってデータをダウンロードしてきました。 f:id:kumawave:20210818210043p:plain f:id:kumawave:20210818210109p:plain

データはこんな感じで説明もあります。

0  39.10205 139.00005   672
0  39.09754 139.00021   736
0  39.09303 139.00038   704
0  39.08852 139.00055   676
0  39.08401 139.00070   669
0  39.07949 139.00087   662
0  39.07498 139.00104   661

f:id:kumawave:20210818210247p:plain

点のデータとしてQGISで直接読み込んでみると

f:id:kumawave:20210818210523p:plain f:id:kumawave:20210818210533p:plain

それらしい位置に表示されますね。

これを扱いやすいようにGeoTiffに変換しようと検索したところ、結局GDALという地理関連の処理ライブラリで変換するのが主な方法のようです。

このYoutubeビデオ

https://www.youtube.com/watch?v=zLNLG0j13Cw&ab_channel=MakingSenseRemotely

を見ながらPythonスクリプトを作成し

#%%
from osgeo import gdal
import pandas as pd
import numpy as np
import os

filepath="jodc-depth500mesh-20210728053445.txt"
df = pd.read_csv(filepath, sep="\s+", header=None)
# Keep only xyz
df.drop(columns=[0], inplace=True)
df.columns=['x','y','value']

df = df.sort_values(by=['y','x'], ascending=[False, True])
df.to_csv("dem.xyz", index=False, header=None, sep=" ",float_format='%.8f')
#demn = gdal.Translate('dem.tif', 'dem.xyz', outputSRS='EPSG:32719')
df.to_csv("dem.csv", index=False)

if os.path.exists('dem.vrt'):
    os.remove('dem.vrt')
    
with open("dem.vrt", 'w') as f:
    f.write("<OGRVRTDataSource>\n \
    <OGRVRTLayer name=\"dem\">\n \
        <SrcDataSource>dem.csv</SrcDataSource>\n \
        <GeometryType>wkbPoint</GeometryType>\n \
        <GeometryField encoding=\"PointFromColumns\" x=\"y\" y=\"x\" z=\"value\"/>\n \
    </OGRVRTLayer>\n \
</OGRVRTDataSource>")

g = gdal.Grid("dem.tif", "dem.vrt", outputSRS="EPSG:4326")
g = None

作成されたTiffファイルをQGISで読み込んだところ

f:id:kumawave:20210818211936p:plain

それらしい位置に表示されました。補間の影響と思われる人工的なパターンが見られますがおそらく補間の方法とパラメータの設定で消えそうな感じです。

スクリプトを作成したというより、Youtubeのビデオにあるスクリプトをコピーしただけですね。作者に感謝です。また、こちらの記事も参考になりました。 500mメッシュ水深データをGISで利用する方法 - 自然環境保全のための周辺技術