福野泰介の一日一創 - create every day

データで遊ぼう!

山口県工業技術センターにて開催「衛星データ解析技術研究会」にてオープンデータ特集。
オープンデータと宇宙データを組み合わせると、ますます楽しく遊べます。

自治体5つ星オープンデータ化を支援するodpなど、5つ星オープンデータの広まりを心待ちにしつつも、いますぐプログラミングして遊べる楽しいこと、はじめましょう!


Tellusで使える2010年3月〜2018年7月末までの降雨量データ、緯度経度共に0.1度の精度で降雨量が分かります。これを使って宇部市、常磐公園や宇部高専あたりの晴天率を調べてみました。

2018年6月の晴天率は、60%とでました。
Tellus OSでもチェックできますが、プログラムにするといろんな場所を自動的にチェックできて便利です。

こちらが計算するのに使った、Jupyter Notebookで動かすPythonプログラム。

TOKEN = "xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx" import os, json, requests, math from skimage import io from io import BytesIO import calendar import matplotlib.pyplot as plt %matplotlib inline def get_image_GSMaP(date, xtile, ytile): sname = "GSMaP" zoom = 5 url = "https://gisapi.tellusxdp.com/{}/{}/{}/{}/{}.png".format(sname, date, zoom, xtile, ytile) print(url) headers = { "content-type": "application/json", "Authorization": "Bearer " + TOKEN } r = requests.get(url, headers=headers) return io.imread(BytesIO(r.content)) def ll2tile(lat, lng, z): w = math.pow(2, z) / 2 yrad = math.log(math.tan(math.pi * (90 + lat) / 360)) return [ (lng / 180 + 1) * w, (1 - yrad / math.pi) * w ] def get_pixel_GSMaP(date, lat, lng): tile = ll2tile(lat, lng, 5) img = get_image_GSMaP(date, math.floor(tile[0]), math.floor(tile[1])) return img[math.floor(tile[0] * 256 % 256)][math.floor(tile[1] * 256 % 256)] def calc_fine_ratio(y, m, lat, lng): cnt = 0 _, ndays = calendar.monthrange(y, m) for d in range(1, ndays + 1): date = "{:0=4}-{:0=2}-{:0=2}".format(y, m, d) rgba = get_pixel_GSMaP(date, lat, lng) #n = (rgba[0] << 24) | (rgba[1] << 16) | (rgba[2] << 8) | rgba[3] n = rgba[3] print(date, n) if n < 255: cnt = cnt + 1 return cnt / ndays calc_fine_ratio(2018, 6, 33.931245,131.277600)

緯度経度からタイル座標への返還は、こちらJavaScriptのプログラムを参考にしています。


PCN山口代表のコヤマさんにも会えました!宇部高専生によるPCN宇部との連携、はじまります!


さくらインターネット、小笠原さん、牟田さんとも楽しいパネルディスカッション。


講演資料オープンデータ「オープンデータと宇宙データで地域活性化(PDF)


日本の衛星、だいち2号(ALOS-2)の1/3スケール模型。この子が宇宙からデータを届けてくれているんですね!
1/1スケール模型、Oculus QuestでぜひVR表示させたいところ!3Dオープンデータ探します。


山口名物、瓦そば!
瓦そばオープンデータがあったら、衛星データと組み合わせて、どんなアプリを作りますか?
こどももおとなも、レッツ、プログラミング!

宇宙データプラットフォームを使った分析はじめ、続いては、お気に入りデバイス Oculus Questと組み合わせます。
2.5mの分解能で宇宙から標高を計測する、ALOS-PRISM。Tellusの開発環境、Jupyther Notebookを使って富士山付近のデータをエクスポートして、VR表示しました。


VR fujisan
マイクラでリアルなマップを再現したい夢、データとプログラミングを使えば簡単に叶います!

作り方
1. Tellusに、Jupyther Notebook か Jupyter Lab を利用申請。
2. 富士山の位置のメッシュを調べる。(タイル座標確認ツール
3. Jupyter Notebook上、Pythonでプログラムをかく

TOKEN = "xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx" import os, json, requests, math from skimage import io from io import BytesIO import matplotlib.pyplot as plt %matplotlib inline def get_image(zoom, xtile, ytile, opacity=1, r=3, g=2, b=1, rdepth=1, gdepth=1, bdepth=1, preset=None): sname = "aw3d30" url = "https://gisapi.tellusxdp.com/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format(sname, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth) if preset is not None: url += '&preset='+preset headers = { "content-type": "application/json", "Authorization": "Bearer " + TOKEN } r = requests.get(url, headers=headers) return io.imread(BytesIO(r.content)) img = get_image(10,906,404) io.imshow(img)

富士山の標高データが画像で表示されます。
4. RGB値が最高になっているようなので、最大値を求めてみる(きっともっとシンプルに書けます)

max = 0 for d in img: n = (d[0] << 16) | (d[1] << 8) | d[2] if n > max: max = n print(max)

結果

12766

富士山の高さは3766mなので、9000ずれた1m単位だと推定。 5. 試しにテキストで表示してみる

for y in range(256): for x in range(256): d = img[x,y] n = (d[0] << 16) | (d[1] << 8) | d[2] n -= 9000 n /= 4000 n = math.floor(n) print(n,end='') print()


それっぽい。
6. JavaScript用に配列データとして出力

h = [] for y in range(256): for x in range(256): d = img[x,y] n = (d[0] << 16) | (d[1] << 8) | d[2] n -= 9000 h.append(n) print(h)

このデータを、fujisan-data.js として保存して、JavaScriptでWebVRを使ったVRアプリをつくって、できあがり!

<!DOCTYPE html> <html> <head> <meta charset="utf-8"/> <title>VR fujisan - data by Tellus</title> <script src="https://aframe.io/releases/0.9.1/aframe.min.js"></script> <script src="https://fukuno.jig.jp/app/webvr/fujisan-data.js"></script> </head> <body> <script>"use strict" window.onload = function() { var scene = document.querySelector('a-scene') const r = 1 const rw = 4 // sampling rw * rw const bw = 256 // org data sampling w and h const w = bw / rw let map = [] for (let i = 0; i < w * w; i++) { const x = i % w const y = (i / w) >> 0 const p = x * rw + y * rw * bw const n = HEIGHT[p] map[i] = n } for (let i = 0; i < map.length; i++) { const x = i % w const y = (i / w) >> 0 const h = map[i] * 0.0005 let box = document.createElement('a-box') box.setAttribute('position', { x: (x - w / 2) * r, y: h / 2 * r - 7, z: (y - w / 2) * r }) box.setAttribute('depth', r) box.setAttribute('width', r) box.setAttribute('height', r * h) scene.appendChild(box) } } </script> <a-scene id="scene"> <a-sky color="#000000"></a-sky> <a-light color="blue" position="1 0.1 0.5"></a-light> <a-light color="cyan" position="0.1 4 2"></a-light> </a-scene> </body> </html>

縮尺は適当なのでリアルな富士山ではないですが、これで怪獣や巨人になった気分で富士山付近を歩けます。
VR x Tellus、楽しいですよ!

links
- VRではじめる現代HTML入門 - Oculus Quest x 福井高専生
- 50行で作るVRマイクラ風 / boxcraft for Oculus Quest with A-Frame
- モバイル時代こそバス! つつじバスロケVRビジュアライズ - リアルタイムオープンデータ x Oculus Quest
- 異次元の分かりやすさ、触って学ぶVR数学、ベジェ曲線編 on Oculus Quest

600万円、これは世界一高い参議院選挙立候補に必要な供託金(市議選は30万円)。
第25回参議院選挙の投票日は7/21。 前回投票日の2016年7月10日の投票所と人の流れをTellusで見てみました。


Tellus ver 1.1 によって追加された、位置情報データ「Profile Passport」を使って、投票日7/10の10時時点(日本時間はUTC+9時間)を品川区の5つ星オープンデータ、投票所と組み合わせて表示した図。 もう一歩分解能がほしいところですが、投票所への人の流れがわかるかも?通常の日曜日との比較なども可能です。


鯖江市の選挙ポスター掲示場、オープンデータと宇宙データ、だいち(ALOS)のセンサーAVRIR-2による植物活性度を重ねたもの。鯖江市の東側、自然豊かな場所にあることが一目瞭然。


Tellusを使った分析は、Jupyter Notebookや、後継サービスJupyter Lab上のPythonで行います。
5つ星オープンデータを集約する opendata.cc のAPI、SPARQLを呼び出し、GeoJSONにしてpostするとTellus上の「取り込みマップ」として表示できます。

API_TOKEN = "xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx" # Tellus マイページ、APIアクセス設定で生成したトークン import json, requests, urllib.parse URL = "https://api.tellusxdp.com/v1/geojson-files" def get(id = ""): return requests.get(URL + "/" + id , headers = { "Authorization": "Bearer " + API_TOKEN }).json() def post(data): headers = { "Authorization": "Bearer " + API_TOKEN, "Content-Type": "application/json" } return requests.post(URL, headers = headers, json = data).json() def delete(id): return requests.delete(URL + "/" + id, headers = { "Authorization": "Bearer " + API_TOKEN }).json() SPARQL_ENDPOINT = "https://sparql.opendata.cc/data/sparql?output=json&query=" def sparql(query): url = SPARQL_ENDPOINT + urllib.parse.quote(query) return requests.get(url).json() def toGeoJSON(res): fs = [] for d in res["results"]["bindings"]: f = { "type": "Feature", "geometry": { "type": "Point", "coordinates": [ d["lng"]["value"], d["lat"]["value"] ] } } if "name" in d: f["properties"] = { "name": d["name"]["value"] } fs.append(f) return { "type": "FeatureCollection", "features": fs } def addLayer(rdftype, layername): gj = toGeoJSON(sparql(""" select ?uri ?name ?lat ?lng { ?uri <http://www.w3.org/1999/02/22-rdf-syntax-ns#type> <""" + rdftype + """>. optional { ?uri <http://www.w3.org/2000/01/rdf-schema#label> ?name. } ?uri <http://www.w3.org/2003/01/geo/wgs84_pos#lat> ?lat. ?uri <http://www.w3.org/2003/01/geo/wgs84_pos#long> ?lng. } """)) #print(gj) post({ "originalName": layername, "transparency": 1, "windowNumber": 1, "data": gj }) addLayer("http://odp.jig.jp/odp/1.0#Polls", "odp 投票所") addLayer("http://odp.jig.jp/odp/1.0#PosterPlace", "odp 選挙ポスター掲示場") print("ok")

rdftypeに設定するデータ種別と対応地区は「5つ星オープンデータ対応一覧表」をご参照ください。


Tellus x オープンデータ、いろいろ遊んでみましょう!

Tweet
クリエイティブ・コモンズ・ライセンス
この作品は「Creative Commons — CC BY 4.0」の下に提供されています。
CC BY 福野泰介 - Taisuke Fukuno / @taisukef / アイコン画像 / プロフィール画像