標高demデータから楕円体高のデータをgdalで作成する方法

なんだか楕円体高のdemデータが欲しくなるときがあります。
そんな時のために取得方法を紹介します。

2022.8.19追追記
Phantom4 RTKの地形認識モード使用時にPPKの場合（D-RTK 2を使用しない）は、高度は気圧計で算出されるため、標高やその日の天気で誤差が大きくなってしまうようです。で、やっぱり楕円体高でDEMを用意する必要があるようです。

↓ジオイド高の分だけ高く飛んだというのは、偶然それぐらいの誤差が出ていただけだった可能性が高いので取り消し

2021.10.19追記
~~Phantom4 RTKの地形認識モードのDSMは、普通に標高で良いのかもしれません....実際に飛ばしたら、ジオイド高の分だけ高く飛んでいました...~~

QGIS 3.20をインストールして最新のgdal（とproj）を使えるようにします。

C:\Program Files\QGIS 3.20.3\share\proj　に入っている jp_gsi_gsigeo2011.tifを作業フォルダにコピーします。

OSGeo4W shellを起動して作業フォルダに移動します。

日本のgeoidモデルのtifをgtx形式に変換します。

gdal_translate -ot Float32 jp_gsi_gsigeo2011.tif jp_gsi_gsigeo2011.gtx

※ジオイドモデルの元データは、国土地理院で提供されています。使用に関する許可申請等は国土地理院のサイトおよびprojファイルに入っているjp_gsi_README.txtを確認ください。

jp_gsi_gsigeo2011.gtxをC:\Program Files\QGIS 3.20.3\share\projのフォルダにコピーします。

標高demデータを緯度経度（epsg:6668）で用意します。

gdalwarpで楕円体高になるようにパラメータを指定して変換します。

gdalwarp -s_srs "+init=epsg:6668 +geoidgrids=jp_gsi_gsigeo2011.gtx" -t_srs epsg:6667 merge.tif merge_ellipsoidal.tif

※ Warning 1: +init=epsg:XXXX syntax is deprecated.　と出ますが気にしません。

こんなことをしなくても gdalwarpで -s_srs をepsg:6697 （標高）-t_srsをepsg:6667（楕円体高）と指定すれば変換できると思ったのですが駄目でした。なんで駄目かを誰か教えてください！

水域などnodataの場所がある場合は気をつけてください！（nodataを0とかにしてるとジオイド高が値に入って、それっぽくなってしまいます。そのままドローンを飛ばすと水域にダイブすることになります。）