オープンデータで航空写真や標高データをPythonで作る方法

Thu, May 14, 2020

${{DZ_TITLE}}$
国土地理より航空写真や地図の画像データがURLとして公開されているのは知っているでしょうか？
地理院タイルと呼ばれていて、タイル状に並んでいる画像のサイトがあります。
任意の範囲の画像を作りには、必要なタイルを取得し貼り合わせ、無駄な部分をトリミングして完成となる。
微妙にめんどくさい作業ですが、Class化しておいたので誰でも簡単に取得できる感じにまとめてみました。

地理院タイルとは

Google mapをほとんどの人は触ったことがあるだろうが、地図モードではなく写真を使った航空写真モードがあるのも大抵の人も知っているでしょう。
日本でも国土地理院が日本の航空写真🗾をタイル画像として公開している。

また、航空写真にとどまらずかなりいろいろな種類のタイル情報があります。（後述）

タイル画像は下の様なURL形式で表されています。
https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png
{z}/{x}/{y}部分を置き換えてみるとURLが完成する。

例えば以下のURLにアクセスし他場合、東京タワーが見えるはずです。
http://cyberjapandata.gsi.go.jp/xyz/ort/17/116415/51623.jpg
2020-05-14-001

画像の使用許諾

もしも本格的にこちらの手法を使う場合、使用許諾等確認してください。
凄く緩いので通常は問題ないですが、場合により引っかかる場合もあるので要注意です。
https://maps.gsi.go.jp/development/ichiran.html

プログラム

import os
import math
import requests
import time
import cv2


class GetImage :
    """
    航空写真をWebデータから生成する

    """

    def long2x_float( long , zoom ):
        """
        Logtitude(経度)より タイル画像の座標を求める。（小数点も活用の必要あり）

        Parameters
        ----------
        long float
            Logtitude(経度)
        zoom int
            拡大値
        Returns
        -------
        x
            タイル画像中の座標(1=タイル1枚分)
        """
        center_point = 2 ** ( zoom + 7 )
        total_pixels = center_point * 2
        pixels_per_long_degree = total_pixels / 360;

        return (center_point + long * pixels_per_long_degree + 0.5) / 256

    def lat2y_float( lat , zoom ):
        """
        Latitude(緯度)より タイル画像の座標を求める。（小数点も活用の必要あり）

        Parameters
        ----------
        lat float
            Latitude(緯度)
        zoom int
            拡大値
        Returns
        -------
        y
            タイル画像中の座標(1=タイル1枚分)
        """
        center_point          = 2 ** ( zoom + 7 )
        total_pixels          = center_point * 2
        pixels_per_lng_radian = total_pixels / (2 * math.pi)
        siny = min(max(math.sin(lat * (math.pi / 180)), -0.9999), 0.9999);

        return       (center_point - 0.5 * math.log((1 + siny) / (1 - siny)) * pixels_per_lng_radian + 0.5) / 256

    def get_xy_float( long , lat , zoom ):
        """
        Logtitude(経度)、Latitude(緯度)より タイル画像の座標を求める。（小数点も活用の必要あり）

        Parameters
        ----------
        long float
            Logtitude(経度)
        lat float
            Latitude(緯度)
        zoom int
            拡大値
        Returns
        -------
        x
            タイル画像中の座標(1=タイル1枚分)
        y
            タイル画像中の座標(1=タイル1枚分)
        """
        return GetImage.long2x_float( long , zoom ) , GetImage.lat2y_float( lat , zoom )

    def make_url_by_xy( x, y, zoom, name ):
        """
        Logtitude(経度)、Latitude(緯度)より タイル画像のURLを取得する

        Parameters
        ----------
        x int
            タイル画像中の座標(1=タイル1枚分)
        y int
            タイル画像中の座標(1=タイル1枚分)
        zoom int
            拡大値
        name str
            イメージ取得URLのAlias名
        Returns
        -------
        image_url
            イメージURL
        """
        dic = {}
        # 必要に応じて追加してください。
        # 下記URLに大量に情報がありますよ。
        # https://maps.gsi.go.jp/development/ichiran.html
        dic['cyberjapandata.airphoto'] ='http://cyberjapandata.gsi.go.jp/xyz/airphoto/{z}/{x}/{y}.png'
        dic['cyberjapandata.ort'] = 'http://cyberjapandata.gsi.go.jp/xyz/ort/{z}/{x}/{y}.jpg'
        dic['cyberjapandata.relief'] = 'https://cyberjapandata.gsi.go.jp/xyz/relief/{z}/{x}/{y}.png'
        # https://maps.gsi.go.jp/development/elevation.html 標高情報
        # 15
        dic['cyberjapandata.dem5a_png'] = "https://cyberjapandata.gsi.go.jp/xyz/dem5a_png/{z}/{x}/{y}.png"
        # 15
        dic['cyberjapandata.dem5b_png'] = "https://cyberjapandata.gsi.go.jp/xyz/dem5b_png/{z}/{x}/{y}.png"
        # 14
        dic['cyberjapandata.dem_png'] = "https://cyberjapandata.gsi.go.jp/xyz/dem_png/{z}/{x}/{y}.png"
        return dic[name].replace('{z}',str(zoom)).replace('{x}',str(x)).replace('{y}',str(y))

    def get_filepath( x , y , zoom , name , dirname ,ext):
        """
        Logtitude(経度)、Latitude(緯度)より キャッシュファイルの保存名を取得
        フォルダ作成も同時に行う

        Parameters
        ----------
        x int
            タイル画像中の座標(1=タイル1枚分)
        y int
            タイル画像中の座標(1=タイル1枚分)
        zoom int
            拡大値
        name str
            イメージ取得URLのAlias名
        Returns
        -------
        image_url
            イメージURL
        """
        cache_dirname = os.path.join( dirname , 'Cache-GetImage', name , str(zoom))
        os.makedirs(cache_dirname, exist_ok=True)

        filename = '%d-%d%s' % ( x , y ,ext)
        filepath = os.path.join( cache_dirname , filename )
        print(filepath)
        return filepath

    def get_file( x , y , zoom , dirname, name, ext , sleep_time = 5) :
        """
        ファイルをダウンロードしてキャッシュに保存する

        Parameters
        ----------
        x int
            タイル画像中の座標(1=タイル1枚分)
        y int
            タイル画像中の座標(1=タイル1枚分)
        zoom int
            拡大値
        dirname str
            保存フォルダ
        name str
            定義名
        sleep_time int
            スリープ時間（サーバーに負荷をかけすぎないように）
        Returns
        -------
        filepath
            取得したファイル名
        """
        filepath = GetImage.get_filepath( x , y , zoom , name, dirname, ext)
        if os.path.exists( filepath ) :
            return filepath

        time.sleep(sleep_time)

        url = GetImage.make_url_by_xy( x , y , zoom ,name)
        print(url)
        res = requests.get( url )
        with open(filepath,'wb+') as fp:
            fp.write(res.content)
        return filepath

    def make_file_by_xy( pos1, pos2, zoom, dirname,  name, ext, out_filepath = None, sleep_time = 5) :        
        """
        画像ファイルを、タイル座標X,Yを元に生成。内部でダウンロードも実施

        Parameters
        ----------
        pos1 list
            タイル座標 X,Y
        pos2 list
            タイル座標 X,Y
        zoom int
            拡大値
        dirname str
            保存フォルダ
        out_filepath str
            出力ファイル名
        name str
            定義名
        sleep_time int
            スリープ時間（サーバーに負荷をかけすぎないように）
        Returns
        -------
        filepath
            取得したファイル名
        """
        x1 = min(pos1[0],pos2[0])
        x2 = max(pos1[0],pos2[0])
        y1 = min(pos1[1],pos2[1])
        y2 = max(pos1[1],pos2[1])
        x1 = int(x1)
        x2 = int(x2)
        y1 = int(y1)
        y2 = int(y2)

        img_h = None
        for x in range( x1, x2 +1 ):
            img_v = None
            for y in range( y1, y2 +1 ):
                filepath = GetImage.get_file( x , y , zoom , dirname, name, ext, sleep_time = 5)
                im1 = cv2.imread( filepath )
                if img_v is None:
                    img_v = im1
                else:
                    img_v = cv2.vconcat([img_v, im1])
            if img_h is None :
                img_h = img_v
            else:
                img_h = cv2.hconcat([img_h, img_v])

        if out_filepath is not None:
            cv2.imwrite(out_filepath, img_h)
        return img_h

    def make_file( pos1, pos2, zoom, out_filepath, name = 'cyberjapandata.ort', ext ='.jpg', max_size = (20,20) , sleep_time = 5) :
        """
        画像ファイルを、Logtitude(経度)、Latitude(緯度)を元に生成。内部でダウンロードも実施

        Parameters
        ----------
        pos1 list
            Logtitude(経度)、Latitude(緯度)
        pos2 list
            Logtitude(経度)、Latitude(緯度)
        zoom int
            拡大値
        out_filepath str
            保存ファイル名
        max_size list
            X , Y の最大タイル数
        sleep_time int
            スリープ時間（サーバーに負荷をかけすぎないように）
        Returns
        -------
        """
        dirname , _ = os.path.split(out_filepath)
        long1 = min(pos1[0],pos2[0])
        long2 = max(pos1[0],pos2[0])
        lat1  = min(pos1[1],pos2[1])
        lat2  = max(pos1[1],pos2[1])
        x1,y1 = GetImage.get_xy_float( long1, lat1, zoom)
        x1_offset = int( ( x1 - int(x1) ) * 256 )
        y2_offset = int( ( y1 - int(y1) ) * 256 )
        x2,y2 = GetImage.get_xy_float( long2, lat2, zoom)
        x2_offset = int( ( x2 - int(x2) ) * 256 )
        y1_offset = int( ( y2 - int(y2) ) * 256 )

        if x2-x1 > max_size[0] :
            raise Exception( "It's over maxsize (lat) : " , x2-x1 , max_size[0])
        if y1-y2 > max_size[1] :
            raise Exception( "It's over maxsize (lng) : " , y1-y2 , max_size[1])

        im = GetImage.make_file_by_xy( (x1,y1), (x2,y2), zoom, dirname, name, ext)
        print(im.shape,x1_offset,y1_offset,x2_offset,y2_offset)

        y2_offset = im.shape[0] - 256 + y2_offset
        x2_offset = im.shape[1] - 256 + x2_offset
        cv2.imwrite(out_filepath, im[y1_offset:y2_offset,x1_offset:x2_offset])

Class定義が完了したので、出力してみましょう。

pos1         = (129.740683, 32.625183)
pos2         = (129.735351, 32.630297)
zoom         = 17
out_filepath = r'D:\temp\out.jpg'
GetImage.make_file( pos1, pos2, zoom, out_filepath )

# Jupyter Notebook環境下で画像を表示
from IPython.display import Image, display_jpeg
display_jpeg(Image(out_filepath))

    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112771/52955.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112771/52956.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112771/52957.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112771/52958.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112772/52955.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112772/52956.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112772/52957.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112772/52958.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112773/52955.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112773/52956.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112773/52957.jpg
    http://cyberjapandata.gsi.go.jp/xyz/ort/17/112773/52958.jpg
    (1024, 768, 3) 51 237 36 35

2020-05-14-002

無事に軍艦島が表示されましたね。

では！

たとえば、rerifを表示したら、こんな感じに表示されます。
こちらは富士山
2020-05-14-003

本当はこちらを知るきっかけになったサイトがあるのですが、とある事情があるのでこちらでは紹介差し控えます。
何故かというと、こちら G??GLE M?Pのデータも同じ方法で取得しているのですが、G??GLEさんそんな事認めていないはずなので…ちょっと紹介できないなぁ…と言う事感じです。

色々いじってチャレンジしてみてください。

標高データのダウンロード

追加報告です。
標高データもメッシュデータとして提供されているようです。

pos1         = (138.572320, 35.488748)
pos2         = (138.864831, 35.290027)
zoom         = 14
out_filepath = r'D:\temp\out.jpg'
GetImage.make_file( pos1, pos2, zoom, out_filepath , name='cyberjapandata.dem_png', ext='.png' , max_size=(27,27))

# Jupyter Notebook環境下で画像を表示
from IPython.display import Image, display_png
display_jpeg(Image(out_filepath))

こんな感じのPNGファイルになります。
なんじゃこりゃって、画像ですね。
保存には必ずPNG等の画像が劣化しないファイルフォーマットを使う必要があります。
色が少しでも違ってしまうと意味をなさくなってしまうためです。
下記のデータはWeb用に画質JPEGで落としてあるので、使い物にならないです。
2020-05-14-004 これを標高に直すには下記の様に計算する必要があるとの事。

x = 216R + 28G + B
x < 223の場合　h = xu
x = 223の場合　h = NA
x > 223の場合　h = (x-224)u

uは標高分解能（0.01m）を表します。
また、無効値（標高タイル（テキスト形式）の「e」に該当する箇所）は(R, G, B)=(128, 0, 0)です。
参照： https://maps.gsi.go.jp/development/demtile.html

標高用のメッシュは、３つのサイトで管理されZoomが15,15,14の固定です。
また、dem5a_pngとdem5b_pngは5mmメッシュ、dem_pngが10mメッシュです。
dem5a_pngとdem5b_pngにはデータのない箇所が多く、特にdem5b_pngは顕著なように思われます。

# https://maps.gsi.go.jp/development/elevation.html 標高情報
# Zoom:15 5mmメッシュ
dic['cyberjapandata.dem5a_png'] = "https://cyberjapandata.gsi.go.jp/xyz/dem5a_png/{z}/{x}/{y}.png"
# Zoom:15 5mmメッシュ
dic['cyberjapandata.dem5b_png'] = "https://cyberjapandata.gsi.go.jp/xyz/dem5b_png/{z}/{x}/{y}.png"
# Zoom:14 10mmメッシュ
dic['cyberjapandata.dem_png'] = "https://cyberjapandata.gsi.go.jp/xyz/dem_png/{z}/{x}/{y}.png"