地理院地図から画像を取得し、任意の緯度経度に文字を入れる方法

Posted On 2021-08-08

7.7k{icon} {views}

Google Map APIの代替として、無料で使える地理院タイルを使ってみました。Pythonのstaticmapと併用すれば、緯度経度から地図画像に文字を入れられました。地理院地図だけでなく空中写真での利用も可能です。

Google Map APIではできるけど

Google Mapのような地図アプリから任意の領域の画像を取得し、適当な地点に文字を入れるという処理を考えます。今までこういうのはパワポでやっていましたが、手動でやるとそこそこ時間がかかります。地図の画像取得や、緯度経度がわかっている地点（例：駅など）への文字入れを自動化できないかというのが今回のコンセプトです。

Webサイトへの埋め込みではなく、画像として地図を取得するスタイルをStatic Mapというそうです。Google MapでもMap Static APIというのがあります。冒頭の「A quick example」というのは、まさに任意の領域の地図を画像として取得し、指定した地点にピンをおいている例ですね（画像はGoogleのサイトより引用）。

はるか昔はAPIキー（認証）なしでもStatic Mapをある程度取れたそうですが、現在は有料のAPIが必須です。Map Static APIは1月につき「2.00 USD per 1000 Map load」と、Googleのサービスにしてはかなり安いです。ただアプリとして作るわけではなく、単に数十～数百のマップ画像に文字入れだけで、このAPIを使うほどでもないなという感はすごいします。

そこで今回無料の代替案を探ってみました。

Static Mapを取得するPythonライブラリ「staticmap」

Pythonライブラリで「staticmap」というのがあります。

https://github.com/komoot/staticmap

Google Mapや地理院地図は地図タイルという、例えば256×256ピクセルの細切れの地図画像を集めた形式で管理しています。本来はこのタイル画像を結合する処理が必要になりますが、staticmapのライブラリを使うとライブラリ側でよしなに結合して必要な領域を切り出してくれます。指定は緯度経度でOKです。

ズームレベル

ライブラリを使う上で緯度経度の指定でOKと書きましたが、縮尺のパラメーターだけは「ズームレベル」という地図タイル特有の値を使用します。

この値はいわゆる「○万分の1」という縮尺と必ずしも対応しません。広く使われてはいるものの、明確に統一された定義ではないそうです。以下では地理院地図の定義を使います。国土地理院のサイトが詳しいです。

地理院タイルについて

簡単に言うと

地球ほぼ全体（北緯・南緯85.0511度以上は除外）をメルカトル投影で1枚の地図として考える。このときがズームレベル0
ズームレベル0の画像を縦横に2等分したものが、ズームレベル1。このとき地球は4枚の地図からなる。
以下同様にズームレベル+1ごとに、縦横2等分する。地球全体の地図の枚数は「4^ズームレベル」のようにべき乗で増えていく。

画像は国土地理院のサイトより引用。ズームレベルが上がるほど拡大された地図になります。

実際のズームレベルってどう見ればいいの？

地理院地図の場合はURLにズームレベルが表示されます。

つまり、地理院地図のサイトを見ていい感じの縮尺になったところで、そのズームレベルをメモすればいいというわけです。

開発者向けに指定した緯度、経度、ズームレベルでタイルを表示するサービスを用意してくれているので、これを活用するのもあり（国土地理院の開発者向けサービス進んでるー）

https://maps.gsi.go.jp/development/tileCoordCheck.html#5/35.362/138.731

地理院地図を画像として取得する

単に地理院地図を画像として取得します。このとき設定するのは「地図の中心の緯度経度、ズームレベル、出力画像の解像度」です。地図タイルの結合はstaticmapライブラリでやってくれます。

参考：Python + staticmap で OpenStreetMap や地理院地図の画像を取得する
https://qiita.com/niwasawa/items/bb86fc1fdf9aeb31d9f0

pip install staticmap

from staticmap import StaticMap
import matplotlib.pyplot as plt

def simple_map():
    # 地理院地図タイル 横1280x縦720
    map = StaticMap(1280, 720, url_template="https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png")
    # 帯広駅を中心にする
    # ズームレベル12、中心画像は「東経, 北緯」→PILのインスタンスとして返ってくる
    img = map.render(zoom=12, center=[143.202056, 42.918])
    # 表示
    plt.imshow(img)
    plt.show()

if __name__ == "__main__":
    simple_map()

map.renderの結果はPILのインスタンスとして返ってきます。PIL（画像処理ライブラリ）を使っている方なら馴染み深いですね。

緯度経度→地図画像の座標の変換をどうするか

ただ「指定した緯度経度に文字を入れたい」という場合は問題があります。staticmapライブラリは、地図上に線を引く機能や、マーカーを置く機能はありますが、文字を入れる機能はありません（v0.5.5時点）。つまりここはPILで自力で入れる必要があります。

問題となるのは、緯度経度から地図画像の座標をどう計算するかです。結論からいうと、staticmapでこの計算はやっているので、ソースから拝借すればどうにかなります。この計算式は国土地理院のとある論文に書いてありました。

WEB ブラウザ上での表示に適した配信地図データの作成技法とその応用
https://www.gsi.go.jp/common/000076319.pdf

staticmapライブラリもこの計算式を踏襲しています。このあたりに緯度経度から変換している式がありますね。

def _lon_to_x(lon, zoom):
    """
    transform longitude to tile number
    :type lon: float
    :type zoom: int
    :rtype: float
    """
    if not (-180 <= lon <= 180):
        lon = (lon + 180) % 360 - 180

    return ((lon + 180.) / 360) * pow(2, zoom)


def _lat_to_y(lat, zoom):
    """
    transform latitude to tile number
    :type lat: float
    :type zoom: int
    :rtype: float
    """
    if not (-90 <= lat <= 90):
        lat = (lat + 90) % 180 - 90

    return (1 - log(tan(lat * pi / 180) + 1 / cos(lat * pi / 180)) / pi) / 2 * pow(2, zoom)

https://github.com/komoot/staticmap/blob/master/staticmap/staticmap.py#L121-L144

ただし、このソースの計算式は「タイル単位」の座標で、ここでの「1離れている＝隣のタイルにある」を意味します。ピクセルに変換するにはもう1手間いります。

その部分がこちらです。

    def _x_to_px(self, x):
        """
        transform tile number to pixel on image canvas
        :type x: float
        :rtype: float
        """
        px = (x - self.x_center) * self.tile_size + self.width / 2
        return int(round(px))

    def _y_to_px(self, y):
        """
        transform tile number to pixel on image canvas
        :type y: float
        :rtype: float
        """
        px = (y - self.y_center) * self.tile_size + self.height / 2
        return int(round(px))

https://github.com/komoot/staticmap/blob/master/staticmap/staticmap.py#L357-L373

StaticMap内の関数にネストして定義されていました。この部分を取り出して一本にまとめてしまいましょう。

変換式を取り出してみた

def lon_to_pixel(lon, map):
    # 経度→タイル番号
    if not (-180 <= lon <= 180):
        lon = (lon + 180) % 360 - 180
    x = ((lon + 180.) / 360) * pow(2, map.zoom)
    # タイル番号→キャンバスの座標
    pixel = (x - map.x_center) * map.tile_size + map.width / 2
    return int(round(pixel))

def lat_to_pixel(lat, map):
    # 緯度→タイル番号
    if not (-90 <= lat <= 90):
        lat = (lat + 90) % 180 - 90
    y = (1 - log(tan(lat * pi / 180) + 1 / cos(lat * pi / 180)) / pi) / 2 * pow(2, map.zoom)
    # タイル番号→キャンバスの座標
    pixel = (y - map.y_center) * map.tile_size + map.height / 2
    return int(round(pixel))

まとめたものがこちら。第2引数はStaticMapのインスタンスに変えてみました。返り値はPILの画像上のピクセル値になります。

実際に文字をプロットしてみた

from staticmap import StaticMap
import matplotlib.pyplot as plt
from math import log, tan, pi, cos
from PIL import ImageFont, ImageDraw

data = [
    ["御影", 42.942419, 142.935739],
    ["芽室", 42.910022,143.048336],
    ["大成", 42.916056,143.069694],
    ["西帯広", 42.924117,143.125425],
    ["柏林台", 42.928853,143.1645],
    ["帯広", 42.918,143.202056],
    ["札内", 42.912697,143.256247],
    ["幕別", 42.907844,143.359197],
    ["利別", 42.934678,143.424369],
    ["池田", 42.921456,143.453136]
]

def lon_to_pixel(lon, map):
    # 経度→タイル番号
    if not (-180 <= lon <= 180):
        lon = (lon + 180) % 360 - 180
    x = ((lon + 180.) / 360) * pow(2, map.zoom)
    # タイル番号→キャンバスの座標
    pixel = (x - map.x_center) * map.tile_size + map.width / 2
    return int(round(pixel))

def lat_to_pixel(lat, map):
    # 緯度→タイル番号
    if not (-90 <= lat <= 90):
        lat = (lat + 90) % 180 - 90
    y = (1 - log(tan(lat * pi / 180) + 1 / cos(lat * pi / 180)) / pi) / 2 * pow(2, map.zoom)
    # タイル番号→キャンバスの座標
    pixel = (y - map.y_center) * map.tile_size + map.height / 2
    return int(round(pixel))

def draw_txt_example():
    # 地理院地図タイル
    map = StaticMap(1280, 720, url_template="https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png")
    # 帯広を中心にする
    global data
    img = map.render(zoom=12, center=[data[5][2], data[5][1]])

    # 日本語フォント
    font = ImageFont.truetype("C://Windows/Fonts/yugothb.ttc", 32)
    draw = ImageDraw.Draw(img)
    for point in data:
        x = lon_to_pixel(point[2], map)
        y = lat_to_pixel(point[1], map)
        # テキストを点に対して中央揃えするためにテキストのピクセルサイズを取得
        str_w, str_h = draw.textsize(point[0], font=font)
        # 駅名の描画
        if x >= 0 and x < img.width and y >= 0 and y < img.height:
            draw.text((x-str_w//2, y-str_h//2), point[0], (64, 64, 255), font=font)

    plt.imshow(img)
    plt.show()

if __name__ == "__main__":
    draw_txt_example()

ピクセルレベルの座標さえ取得できればあとはただの画像処理なので楽ちんですね。エリア外の駅も除外できていますし、これが無料でできるのは素晴らしいです。

空中写真にも文字入れできる

実は地理院タイルは相当充実していて、地理院地図以外でもできます。StaticMapの初期化の「url_template」を変えるだけです。URL仕様は地理院地図のサイトを見てください。ものによっては別途利用申請が必要なのもありますので、規約は注意してください。

地理院タイル一覧
https://maps.gsi.go.jp/development/ichiran.html

傾斜量図

URL：https://cyberjapandata.gsi.go.jp/xyz/slopemap/{z}/{x}/{y}.png

雑感

空中写真も年代別に機械的に取れるのびっくりしましたし（もしかしたら解像度があまりよくないかもしれない）、標高データは3D的な活用ができそうで面白そう。使い方によってはGoogle Mapより便利なレベルで活用できそうで、楽しそうなサービスでした。

Shikoan's ML Blogの中の人が運営しているサークル「じゅ～しぃ～すくりぷと」の本のご案内

技術書コーナー

北海道の駅巡りコーナー

Tags:Google Map, staticmap, 地理院地図, 空中写真, 空間情報