Mapboxからタイル画像を取得して一枚の画像に合成する

LabsのツールではなくPythonから直接Mapboxから画像をダウンロードする方法。この方法を使うと解像度に限度がなくなる。

環境の準備
コード

環境の準備

Mapboxのアカウント登録して、キーを取得しておく

Python本体のインストール
※インストール時に「Add Python to PATH」をチェックしておく

以下のパッケージをインストールしておく

Pillow：画像処理ライブラリ
Mercantile：地理座標系のライブラリ
Requests：HTTPライブラリ

コマンドプロンプトを開き、以下のコードを順に実行してインストールしていく

pip install pillow
pip install mercantile
pip install requests

コード

mapbox_combine_tile_images.pyのようなテキストファイルを適当に作成し、中身にコードを書く。

コードの流れ

・座標（緯度・経度）と解像度を設定する
・座標から該当するタイル番号を取得する
・必要なタイル数と番号を計算する
・Mapboxサーバーからタイル画像を取得する
・タイル画像を1枚にまとめる

#
# Mapboxからタイル画像を取得して一枚の画像に合成する
#
import math
access_token = 'ここにmapboxのキーを入力する'

# 緯度経度、解像度、タイル数を入力する
lat_lng = [36.17106700090759, 138.17289817315498]
z = 14          # 解像度（1:地球を2(=2^1)分割、2:地球を4(=2^2)分割・・16:地球を65536(=2^16)分割）
size_tile = 16   # タイル数

print('lat: ' + str(lat_lng[0]))
print('lng: ' + str(lat_lng[1]))

# 該当するタイル番号を取得する
import mercantile
center_tiles = mercantile.tile(lat_lng[1], lat_lng[0], z)

# タイルの範囲を設定する
x_tile_range = [int(center_tiles.x - size_tile/2 + 1),int(center_tiles.x + size_tile/2)]
y_tile_range = [int(center_tiles.y - size_tile/2 + 1),int(center_tiles.y + size_tile/2)]

#
# 緯度からタイルの長さを計算する
#
R = 6378.137    # 地球の半径（km）
C = R * 2 * math.pi * 1000  # 赤道の円周
rad = math.pi/2 * (lat_lng[0]/90.0)     # 緯度をラジアンの角度へ
S = C * math.cos(rad) / math.pow(2, z)  # 緯度における円周を地図の解像度で割り、タイルあたりの長さを求める

print('Zoom: ' + str(z))
print('Rezolution: ' + str(size_tile * 512) + 'px x '+ str(size_tile * 512)+'px')
print('width(m): ' + str(S * size_tile) + '(m) x ' 'height(m): ' + str(S * size_tile) + '(m)')
print('tile width(m): ' + str(S) + '(m)')
print('pixel width(m): ' + str(S/512) + '(m)')

# 左の座標を計算する
bounds_west_north = mercantile.bounds(x_tile_range[0], y_tile_range[0], z)
bounds_east_south = mercantile.bounds(x_tile_range[1], y_tile_range[1], z)

print('')
print('west:' + str(bounds_west_north.west))
print('east:' + str(bounds_east_south.east))
print('north:' + str(bounds_west_north.north))
print('south:' + str(bounds_east_south.south))

print('')
print('Tile X Range: ' + str(x_tile_range))
print('Tile Y Range: ' + str(y_tile_range))
print('Tiles: ' + str(x_tile_range[1] - x_tile_range[0] + 1) + ' x ' + str(y_tile_range[1] - y_tile_range[0] + 1))

print("Press Any Key to Continue...")
input()

#
# Mapboxサーバーから衛星画像のタイルを取得する
#
import requests
import shutil
import os

# ディレクトリが存在していなければ作成する
if not os.path.exists("elevation_images"):
    os.makedirs("elevation_images")

if not os.path.exists("satellite_images"):
    os.makedirs("satellite_images")

if not os.path.exists("composite_images"):
    os.makedirs("composite_images")

# テキスト書き出し
file = open('composite_images\\' + 'z' + str(z) + '.txt', 'w')
file.write('Zoom: ' + str(z) + "\n")
file.write('Rezolution: ' + str(size_tile * 512) + 'px x '+ str(size_tile * 512)+'px' + "\n")
file.write('width(m): ' + str(S * size_tile) + '(m) x ' 'height(m): ' + str(S * size_tile) + '(m)' + "\n")
file.write('west:' + str(bounds_west_north.west) + "\n")
file.write('east:' + str(bounds_east_south.east) + "\n")
file.write('north:' + str(bounds_west_north.north) + "\n")
file.write('south:' + str(bounds_east_south.south) + "\n")
file.close

# タイル画像をダウンロードしていく処理
for i,x in enumerate(range(x_tile_range[0],x_tile_range[1]+1)):
    for j,y in enumerate(range(y_tile_range[0],y_tile_range[1]+1)):
        print(x,y)

        # ハイトマップを取得する
        r = requests.get('https://api.mapbox.com/v4/mapbox.terrain-rgb/'+str(z)+'/'+str(x)+'/'+str(y)+'@2x.pngraw?access_token='+access_token, stream=True)
        with open('./elevation_images/' + str(i) + '.' + str(j) + '.png', 'wb') as f:
            r.raw.decode_content = True
            shutil.copyfileobj(r.raw, f)  
        
        # 衛星画像を取得する
        r = requests.get('https://api.mapbox.com/v4/mapbox.satellite/'+str(z)+'/'+str(x)+'/'+str(y)+'@2x.png?access_token='+access_token, stream=True)
        with open('./satellite_images/' + str(i) + '.' + str(j) + '.png', 'wb') as f:
            r.raw.decode_content = True
            shutil.copyfileobj(r.raw, f)

print("Press Any Key to Continue...")
input()

#
# 取得した画像を合成して一枚の大きな画像にする
#
import math
from PIL import Image
from os import listdir
from os.path import isfile, join

for img_name in ['elevation','satellite']:
    # ディレクトリ内の画像ファイルのリストを取得
    image_files = ['./'+img_name+'_images/' + f for f in listdir('./'+img_name+'_images/')]

    images = [Image.open(x) for x in image_files]

    edge_length_x = x_tile_range[1] - x_tile_range[0] + 1
    edge_length_y = y_tile_range[1] - y_tile_range[0] + 1
    edge_length_x = max(1,edge_length_x)
    edge_length_y = max(1,edge_length_y)
    width, height = images[0].size

    total_width = width * edge_length_x
    total_height = height * edge_length_y

    composite = Image.new('RGB', (total_width, total_height))

    y_offset = 0
    for i in range(0,edge_length_x):
        x_offset = 0
        for j in range(0,edge_length_y):
            tmp_img = Image.open('./'+img_name+'_images/' + str(i) + '.' + str(j) + '.png')
            composite.paste(tmp_img, (y_offset,x_offset))
            x_offset += width
            
        y_offset += height

    composite.save('./composite_images/'+img_name+'_'+str(int(S * size_tile))+'m'+'_z'+str(z)+'.png')

print("Press Any Key to Finish...")
input()