Cartopy で OpenStreetMap を使う - ぽぷらの散歩道

　Basemap から Cartopy へ引っ越し -2- で解決したい問題として、もう少し解像度の高い地図が欲しい、と書きました。今回、OpenStreetMap を使ってみました。

１．OpenStreetMap とは

２．tube_stations example - gallery in Cartopy -

３． A non-empty list of tiles should be provided to merge.

４．OpenStreetMap API のリクエストにはちょっとコツがあるらしい

５．気象庁 MSM データを重ねたものの…

１．OpenStreetMap とは

　2004 年にイギリスで始まった、みんなで地図情報を作り、使えるようにするプロジェクトです。アカウントを作ると、自分で地図上に書き込みたい情報を編集できるようにもなるそうです。日本バージョンもあるみたいです。

www.openstreetmap.org

openstreetmap.jp

２．tube_stations example - gallery in Cartopy -

　Cartopy gallery にも、OpenStreetMap を使ってロンドンの Waterloo 駅周辺の地下鉄の駅を描いた素敵なサンプルがあったので、まずこれを眺めます。

　結論として、以下の code が地図のみ読み込むものとなることが分かりました。実行結果は下の通りです（gallery よりも、少しだけ範囲を広げています）。

import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from cartopy.io.img_tiles import OSM

imagery = OSM()
ax = plt.axes(projection=imagery.crs)
ax.add_image(imagery, 14)　　　　　　　　　　　　# 拡大率 14
ax.set_extent((-0.15, -0.09, 51.49, 51.52))    # 緯度経度

plt.show()

f:id:ohigehige:20200923231421p:plain

３． A non-empty list of tiles should be provided to merge.

　不思議なことに、code をちょこちょこいじっているうちに、動いていた code でValueError: A non-empty list of tiles should be provided to merge.
というメッセージが出るようになりました。不思議なことにそのうちまた動くようになったりします。メッセージについて検索してみると、それなりに話題になってもいるものでした。

４．OpenStreetMap API のリクエストにはちょっとコツがあるらしい

　詳しい理由は、少なくとも私には残念ながら分かりませんが、OpenStreetMap を読み込む時、ちょっと工夫をすることでメッセージの発生を防げるようです。色々解決法は議論されていますが、以下のサイトが実例と共に、素晴らしく分かりやすいです。

makersportal.com

import io
from urllib.request import urlopen, Request
from PIL import Image

def image_spoof(self, tile): 
　　　　# this function pretends not to be a Python script
 url = self._image_url(tile) 
　　　　# get the url of the street map API
 req = Request(url) 
　　　　# start request
 req.add_header('User-agent','Anaconda 3') 
　　　　# add user agent to request
 fh = urlopen(req) 
 im_data = io.BytesIO(fh.read()) 
　　　　# get image
 fh.close() 
　　　　# close url
 img = Image.open(im_data) 
　　　　# open image with PIL
 img = img.convert(self.desired_tile_form) 
　　　　# set image format
 return img, self.tileextent(tile), 'lower' 
　　　　# reformat for cartopy

cimgt.OSM.get_image = image_spoof 
　　　　# reformat web request for street map spoofing
osm_img = cimgt.OSM() 
　　　　# spoofed, downloaded street map

　先に挙げたサイトのコードをお借りしていますが、地図を読み込む code に上の code を付けて読み込むことで、Error メッセージの発生はなくなります。

５．気象庁 MSM データを重ねたものの…

　解像度の高い地図が得られることが分かったところで、気象庁 MSM の風データを重ねてみることにしました。使用データは、Basemap から Cartopy へ引っ越し -3- でも使用した、2020.08.24 台風8号（Bavi）06Z（午後3時）の風データです。下の図が、その結果になります。

f:id:ohigehige:20200923234702p:plain

　あれ・・・、何で風のデータが無いんでしょう。正確には、1つしか描かれていません。ちょっと考えて我ながらガックリしました。

　「MSM は 5 km メッシュのデータしかそもそもない！」

のでした。今年の春にオープンした那覇の 2 本目の滑走路データは地図に反映されているものの、この地図の解像度は MSM データには細か過ぎました。当日は、同じ南西の風でも、18L の方が 18R よりも少し横風成分が少なく、着陸できると言っている飛行機が多い時間帯があったようだったので、それがデータで確認できるかと思ったのですが…。

　実際、MSM 風データの細かさは下の図程度だったのです。

f:id:ohigehige:20200923213040p:plain