世界地図の陰影図作成
2014年07月22日
地球全体を対象とした数値標高モデルのダウンロード
数値標高モデルの調査
海域を含む地球全体を対象とした、標高データが必要になり調査した。
目的は社内で作成している世界地図カレンダーに、陰影をつけてみること。
世界地図カレンダーの範囲が、北緯80度+α、南緯60度弱、西経20度+α、東経340度+αで、日本を中心として東西に重複がある。
候補として挙げたのが、以下のDEM。
- ETOPO1
http://www.ngdc.noaa.gov/mgg/global/global.html
無償。出典表示。
約1.8km間隔。北緯90度～南緯90度。
氷表面（南極とグリーンランド）と岩盤面の２種類。netCDF、binary、xyz、georeferenced tiff形式。
１ファイルで全域。 - SRTM30Plus
http://topex.ucsd.edu/WWW_html/srtm30_plus.html
の中段にあるFTP SRTM30_PLUS from Sandwell Lab digital grided data
ftp://topex.ucsd.edu/pub/srtm30_plus/
無償。商業目的は問合せ。出典表示。
陸域は、SRTM30と高緯度地域にGTOPO30を使用。約1km間隔。
海域は、NGDCのCoastal Relief Modelなどを使用。約1.8km間隔。北緯81度～南緯81度。 - SRTM
http://www2.jpl.nasa.gov/srtm/
の中段にあるSRTM V2 released
http://dds.cr.usgs.gov/srtm/
無償。出典表示。
SRTM3：約90m間隔、高さ精度10m。北緯60度～南緯56度。
SRTM30：約900m間隔、北緯90度～南緯60度。
海域が含まれない。目的の用途には精度が良すぎる。 - ASTER GDEM
https://www.jspacesystems.or.jp/library/archives/ersdac/GDEM/J/index.html
の上段にあるASTER GDEMのロゴ
http://gdem.ersdac.jspacesystems.or.jp/
無償。登録制。出典表示。約30m間隔、高さ精度7m～14m。北緯83度～南緯83度。GeoTIFF形式。
海域が含まれない。目的の用途には精度が良すぎる。 - WorldDEM
http://www.astrium-geo.com/jp/5725-worlddem
有料。約12m間隔、高さ精度2m(相対)、4m(絶対)。
海域が含まれない。目的の用途には精度が良すぎる。 - 全世界デジタル3D地形データ
http://alos-world3d.jp/
有料。約5m間隔、高さ精度5m。
海域が含まれない。目的の用途には精度が良すぎる。
数値標高モデルの選定
海域を含む条件だけで、２つに絞り込まれ、１ファイルで全域をカバーする取り扱いやすさでETOPO1を利用する事とした。
精度も今回の目的には十分。
氷表面と岩盤面はとりえず岩盤面を選択。
形式は取り扱いやすいGeoTIFF形式を選択し、ダウンロード。
ArcGIS users - use binary and convert using 'Float to Raster'と記してあったが、変換の手間を省くため。
数値標高モデルの加工
数値標高モデルの確認
ダウンロードしたGeoTIFFファイル（ETOPO1_Bed_g_geotiff.tif）をArcGIS10.0で開いて気が付いた。
当然だが、日本が中心になっていない。
他の形式や、ArcGIS10.2やSIS7.1、QGIS2.4で開いても変わるかけが無い。
GDALINFOで情報確認。
gdalinfo D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff.tif
Size is 21601, 10801
Corner Coordinates:
Upper Left (-180.0083333, 90.0083333)
Lower Left (-180.0083333, -90.0083333)
Upper Right ( 180.0083333, 90.0083333)
Lower Right ( 180.0083333, -90.0083333)
Center ( 0.0000000, 0.0000000)
http://www.ngdc.noaa.gov/mgg/global/global.html
の上段、左側にあるExtract Custom Grid
http://maps.ngdc.noaa.gov/viewers/wcs-client/
上記から、必要な範囲をドラッグで指定、座標を入力して指定して試すが、選択範囲が広いとメッセージが出てきた。
試しに日本周辺に範囲を狭めて見るとダウンロードリンクが出て、ダウンロード出来たが、これでは意味が無い。
最初にダウンロードしたGeoTIFFファイル（ETOPO1_Bed_g_geotiff.tif）を加工することにした。
ETOPO1を横並びに２つマージ
最初にダウンロードしたGeoTIFFファイル（ETOPO1_Bed_g_geotiff.tif）を横並びに２つマージさせる。
日本が中心で、東西に重複がある陰影図を作成するため。
東側に１つ分移動させ配置した時の座標を、GDALコマンドで与えてみた。360度を超える設定となり上手くいくか試した。
gdal_translate -of GTiff ^
-gcp 0 0 180.0083333 90.0083333 ^
-gcp 21601 0 540.0083333 90.0083333 ^
-gcp 0 10801 180.0083333 -90.0083333 ^
-gcp 21601 10801 540.0083333 -90.0083333 ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_gcp.tif
測地系を定義。オリジナルのETOPO1と同じWGS84。WGS84のEPSGコードは4326。
gdalwarp -t_srs EPSG:4326 ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_gcp.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_warp.tif
GDALINFOで情報確認。小数点以下にズレが出ているが、後から切り抜くのでOKとした。
gdalinfo D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_warp.tif
Size is 21601, 10801
Corner Coordinates:
Upper Left ( 180.0020833, 90.0083333) (180d 0' 7.50"E, 90d 0'30.00"N)
Lower Left ( 180.0020833, -90.0016665) (180d 0' 7.50"E, 90d 0' 6.00"S)
Upper Right ( 540.005, 90.008) (Invalid angle, 90d 0'30.00"N)
Lower Right ( 540.005, -90.002) (Invalid angle, 90d 0' 6.00"S)
Center ( 360.004, 0.003) (360d 0'13.50"E, 0d 0'12.00"N)
横並びに２つのGeoTIFFファイルをマージ。
gdal_merge -of GTiff ^
-o D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_merge.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_warp.tif
世界地図カレンダーの範囲に切り抜き
切り抜く範囲のポリゴンシェープファイルを作成する。
使い慣れたArcGISでポリゴンシェープファイルを新規作成。
- カタログで右クリック-新規作成-シェープファイル-ファイル名、ポリゴン、座標系を設定。
- エディタ-編集の開始。
- まずは適当にポリゴンを作成（作図ツールの四角形など利用）。
- スケッチプロパティを開いて、四隅の座標値を変更。
- 左上 -30.0 85.0
- 右上 350.0 85.0
- 右下 350.0 -70.0
- 左下 -30.0 -70.0
- スケッチ終了。
- エディタ-編集の保存、編集の終了。
上記で作成したポリゴンシェープファイルの範囲で、マージしたGeoTIFFファイルを切り抜き。
gdalwarp -q -cutline ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_E30W350N85S70.shp ^
-crop_to_cutline -of GTiff ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_merge.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_clip.tif
GDALINFOで情報確認。
gdalinfo D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_clip.tif
Size is 2280, 9300
Corner Coordinates:
Upper Left ( -30.0000000, 85.0000000) ( 30d 0' 0.00"W, 85d 0' 0.00"N)
Lower Left ( -30.0000000, -70.0000000) ( 30d 0' 0.00"W, 70d 0' 0.00"S)
Upper Right ( 350.000, 85.000) (350d 0' 0.00"E, 85d 0' 0.00"N)
Lower Right ( 350.000, -70.000) (350d 0' 0.00"E, 70d 0' 0.00"S)
Center ( 160.0000000, 7.5000000) (160d 0' 0.00"E, 7d30' 0.00"N)
上記の座標範囲のままでは、ミラー図法へ投影変換した時に、不具合となる。
そのため、東西の範囲を-180度から180度の360度範囲に変更する。
gdal_translate -of GTiff ^
-gcp 0 0 -180.0000000 85.0000000 ^
-gcp 22800 0 180.0000000 85.0000000 ^
-gcp 0 9300 -180.0000000 -70.000000000 ^
-gcp 22800 9300 180.0000000 -70.0000000 ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_clip.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_clip_gcp.tif
測地系を定義。オリジナルのETOPO1と同じWGS84。WGS84のEPSGコードは4326。
gdalwarp -t_srs EPSG:4326 ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_gcp.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_warp.tif
ミラー図法へ投影変換
世界地図カレンダーの投影がミラー図法であるため、GeoTIFFファイルをWGS84からミラー図法へ投影変換する。
ミラー図法のEPSGコードは54003。
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:54003 ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_clip_warp.tif ^
D:\Home\DEM\ETOPO1\Bedrock\grid\ETOPO1_Bed_g_geotiff_clip_warp_miller.tif
陰影効果の付与
QGIS2.4のラスタ-解析-DEM（地形モデル）の、陰影図モードで陰影効果を付与。
GDALコマンドでは以下のようになる。
gdaldem hillshade ^
D:/Home/DEM/ETOPO1/Bedrock/grid/ETOPO1_Bed_g_geotiff_clip_warp_miller.tif ^
D:/Home/DEM/ETOPO1/Bedrock/grid/ETOPO1_Bed_g_geotiff_clip_warp_miller_hiishade_z10.tif ^
-z 10.0 -s 1.0 -az 315.0 -alt 45.0 -of GTiff