国土地理院が提供する標高タイルをダウンロードしてGeoTiff形式にする簡単なPythonプログラムをメモ用に置いておく。

標高タイルとは

https://maps.gsi.go.jp/development/demtile.html

ここで解説されている、国土地理院が提供するウェブメルカトルタイルのそれぞれの点の標高値が取得できるヤツ。基盤地図情報などとまとめて提供される数値標高モデルとの関係（点の密度はどっちが高いとかね、）はまた別ページで考察したいと思うが、アカウント登録も不要だしはっきり言ってコッチの方が使いやすい。

Pythonコード

基本的にはタイルの標高値でピクセルを埋めるが優先度はDEM5A>DEM5Bとなっている。DEM5Aはレーザ測量で街区や緩斜エリアに多く、DEM5Bは写真測量で山間部に多い。解説はしない。

#!python3
# -*- coding: utf-8 -*-

import os
import sys
import codecs
import csv
import requests
import math
import numpy as np
import osgeo.gdal,osgeo.osr

def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2**zoom
  xtile = int((lon_deg + 180.0)/360.0*n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1/math.cos(lat_rad)))/math.pi)/2.0*n)
  return xtile,ytile

def num2deg(xtile, ytile, zoom):
  n = 2**zoom
  lon_deg = xtile/n*360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi*(1-2*ytile/n)))
  lat_deg = math.degrees(lat_rad)
  return lat_deg,lon_deg

def main(argv):
  try:
    output_filename = argv[1]
    lat = float(argv[2])
    lon = float(argv[3])
    hsize = int(argv[4])
    vsize = int(argv[5])
  except IndexError:
    print('usage: python %s [output_filename] [latitude] [longitude] [hsize] [vsize]' % (argv[0]))
    quit()

  nodata = -1

  z = 15
  x,y = deg2num(lat, lon, z)

  raster = np.zeros((vsize*256, hsize*256), dtype=np.float64)

  headers = {'User-Agent' : 'Mozilla/5.0 (Windows NT 10.0; WOW64; rv:54.0) Gecko/20100101 Firefox/54.0'}

  for i in range(vsize):
    for j in range(hsize):
      tx = x - hsize//2 + j
      ty = y - vsize//2 + i

      url5a = 'http://cyberjapandata.gsi.go.jp/xyz/dem5a/%d/%d/%d.txt' % (z,tx,ty)
      url5b = 'http://cyberjapandata.gsi.go.jp/xyz/dem5b/%d/%d/%d.txt' % (z,tx,ty)
      print(i,j,url5a)
      print(i,j,url5b)
      response1 = requests.get(url5a, headers=headers)
      response2 = requests.get(url5b, headers=headers)

      dem5a = list(csv.reader(response1.text.splitlines()))
      dem5b = list(csv.reader(response2.text.splitlines()))

      if len(dem5a) == 256:
        for p in range(256):
          for q in range(256):
            if dem5a[p][q] == 'e' : raster[i*256 + p][j*256 + q] = nodata
            else: raster[i*256 + p][j*256 + q] = float(dem5a[p][q])
      else:
        print('error dem5a',i,j)
        for p in range(256):
          for q in range(256):
            raster[i*256 + p][j*256 + q] = nodata

      if len(dem5b) == 256:
        for p in range(256):
          for q in range(256):
            if raster[i*256 + p][j*256 + q] == nodata:
              if not dem5b[p][q] == 'e': raster[i*256 + p][j*256 + q] = float(dem5b[p][q])
      else: print('error dem5b',i,j)

  lat0,lon0 = num2deg(x-hsize//2,y-vsize//2,z)
  lat1,lon1 = num2deg(x+hsize//2+1,y+vsize//2+1,z)

  print(lat0,lon0)
  print(lat1,lon1)

  trans = [lon0, (lon1-lon0)/(hsize*256), 0, lat0, 0, (lat1-lat0)/(vsize*256)]

  srs = osgeo.osr.SpatialReference()
  srs.ImportFromEPSG(4326)

  driver = osgeo.gdal.GetDriverByName('GTiff')
  output = driver.Create(output_filename, hsize*256, vsize*256, 1, osgeo.gdal.GDT_Float32)
  output.GetRasterBand(1).WriteArray(raster)
  output.GetRasterBand(1).SetNoDataValue(nodata)
  output.SetGeoTransform(trans)
  output.SetProjection(srs.ExportToWkt())
  output.FlushCache()

if __name__ == '__main__':
  main(sys.argv)

免責事項とライセンス

プログラムは各自の責任のもとで利用されるものとし、何が起こっても一切の責任は負えません。実験や趣味に用いる分には自由に参考にしてもらって構いませんが、何かしらの業務やプロダクト生成の過程でこのプログラムを利用する場合はContact Informationから著作者に通知してください。