Pythonで地理院地図を画像としてファイルに保存

2020-02-07 23:292020-02-07 23:33

地理院地図はウェブページから普通に閲覧は出来るが、ウェブマップサービスとしてタイル座標からそれぞれの画像にアクセスすることもできる。これを使ってPythonから画像を取得して1枚の画像にするプログラムを置いておく。

Pythonコード

ワールドファイルは生成しないのでただの画像としてしか使えない。QGISからxyzタイル使って保存したほうが画面見ながら範囲指定して保存できるしワールドファイル作れるし再投影できるしこっちはオワコン。

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

import sys
import math
import csv
import os.path
import requests
from io import BytesIO
from PIL import Image

def tile2latlon(x, y, z):
    lon = (x / 2.0**z) * 360 - 180
    mapy = (y / 2.0**z) * 2 * pi - pi
    lat = 2 * math.atan(e ** (- mapy)) * 180 / math.pi - 90
    return [lon,lat]

def latlon2tile(lon, lat, z):
    x = int((lon / 180 + 1) * 2**z / 2)
    y = int(((-math.log(math.tan((45 + lat / 2) * math.pi / 180)) + math.pi) * 2**z / (2 * math.pi)))
    return [y,x]

def execute(output, type, zoom, latitude, longitude, hsize, vsize):
  session = requests.session()
  headers = {'User-Agent' : 'Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:55.0) Gecko/20100101 Firefox/55.0'}
  basename, ext = os.path.splitext(output)
  y0, x0 = latlon2tile(longitude, latitude, zoom)
  width = 2*hsize + 1
  height = 2*vsize + 1
  image = Image.new("RGB", (256*width,256*height), (255,255,255))
  for i in range(width):
    for j in range(height):
      x = x0 - hsize + i
      y = y0 - vsize + j
      print('downloading tile at '+str([x,y]))
      response = session.get('http://cyberjapandata.gsi.go.jp/xyz/%s/%d/%d/%d%s' % (type, zoom, x, y, ext), headers=headers)
      subimage = Image.open(BytesIO(response.content))
      image.paste(subimage, (256*i,256*j))
  image.save(output)

def main(argv):
  if len(argv) < 8:
    print('usage: python %s [output] [type] [z] [latitude] [longitude] [hsize] [vsize]' % argv[0])
  else:
    try:
      output = argv[1]
      type = argv[2]
      zoom = int(argv[3])
      latitude = float(argv[4])
      longitude = float(argv[5])
      hsize = int(argv[6])
      vsize = int(argv[7])
    except ValueError:
      print('argument error')
      quit()
    execute(output, type, zoom, latitude, longitude, hsize, vsize)

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