カテゴリー: Tellus

Tellus ASNARO-1

Tellus

https://www.tellusxdp.com/ja/

Result

Original data provided by NEC
Original data provided by NEC

Code

import os, json, requests, math
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
TOKEN = "Input your token"

def get_ASNARO_scene(min_lat, min_lon, max_lat, max_lon):
	url = "https://gisapi.tellusxdp.com/api/v1/asnaro1/scene"     	+ "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
	headers = {
    	"content-type": "application/json",
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return r.json()

scenes = 0
# Orig
# scenes = get_ASNARO_scene(20.425278, 122.933611, 45.557222, 153.986389)
# Fuji
scenes = get_ASNARO_scene(35.365, 138.705, 35.36, 138.71)
# Tokyo Tower
# scenes = get_ASNARO_scene(35.6505805,139.7402442, 35.6585805,139.74324421)
# Sky Tree　35.710067,139.8085064
# scenes = get_ASNARO_scene(35.7100627,139.8085117, 35.7200627,139.8085117)
# Tokyo Station　35.6812405,139.7649308
# scenes = get_ASNARO_scene(35.6812405,139.7649308, 35.6812405,139.7659308)
# home 35.6034393,139.320697
# scenes = get_ASNARO_scene(35.6034393,139.315697, 35.7234393,139.39697)

print(len(scenes))
print(scenes[0]['thumbs_url'])
      
ext_scene = scenes[0]
img_thumbs = io.imread(ext_scene['thumbs_url'])
print(len(img_thumbs))
io.imshow(img_thumbs)
def get_ASNARO_scene(min_lat, min_lon, max_lat, max_lon):
	url = "https://gisapi.tellusxdp.com/api/v1/asnaro1/scene"     	+ "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
	headers = {
    	"content-type": "application/json",
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return r.json()

def get_tile_num(lat_deg, lon_deg, zoom):
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

def get_ASNARO_image(scene_id, zoom, xtile, ytile):
    url = " https://gisapi.tellusxdp.com/ASNARO-1/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

def get_NxN_series_image(scene_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """切り出したタイル画像を結合する"""
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            row.append(get_ASNARO_image(scene_id, zoom, topleft_x + x, topleft_y + y))
        img.append(np.hstack(row))
    return  np.vstack(img)

def get_4x4(lat, log, zoom):
    xtile, ytile = get_tile_num(lat, log, zoom)
    img_4x4 = get_NxN_series_image(scenes[0]["entityId"], zoom, xtile, ytile, 4, 4)
    return img_4x4 

# 東京タワー 35.6585805,139.7432442
# xtile = 232830
# ytile = 103246
zoom=18
lat = 35.6585805  
log = 139.7432442 

scenes = get_ASNARO_scene(lat, log, lat + 0.1, log + 0.1)

(xtile, ytile) = get_tile_num(scenes[0]['clat'], scenes[0]['clon'], zoom)

print(zoom)
print(xtile, ytile)

# img = get_ASNARO_image(scenes[0]['entityId'], zoom, xtile, ytile)
# io.imshow(img)

img_4x4 = get_4x4(lat, log, zoom)

# io.imshow(img_4x4)

import cv2

def gamma(_img):
    # テーブルを作成する。
    table = np.clip(np.arange(256)*2.5,0, 255)
    # [0, 255] でクリップし、uint8 型にする。
    table = np.clip(table, 0, 255).astype(np.uint8)
    return cv2.LUT(_img, table)


from PIL import Image
dst = gamma(img_4x4)

pil_img = Image.fromarray(dst)

pil_img.save('tt.png')

from IPython.display import Image
Image(url="tt.png")

Tellus SLATS

Tellus

https://www.tellusxdp.com/ja/

Result

Original data provided by JAXA
Original data provided by JAXA

Code

import os, json, requests, math
import numpy as np
import dateutil.parser
from datetime import datetime
from datetime import timezone
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt

TOKEN = "Input your token"

def get_tsubame_scene(min_lat, min_lon, max_lat, max_lon):
    url = "https://gisapi.tellusxdp.com/api/v1/tsubame/scene" \
        + "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return r.json()

tsubame_scenes = get_tsubame_scene(35.524375, 139.585347, 35.724375, 140.005347)

print(len(tsubame_scenes))
print(tsubame_scenes[0])

ext_scene = tsubame_scenes[0]

# 表示
io.imshow(io.imread(tsubame_scenes[0]["thumbs_url"]))
from PIL import Image

def get_tile_num(lat_deg, lon_deg, zoom):
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

def get_tsubame_image(scene_id, zoom, xtile, ytile):
    url = " https://gisapi.tellusxdp.com/tsubame/{}/{}/{}/{}.png".format(scene_id, zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

def get_slat_series_image(scene_id, zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """切り出したタイル画像を結合する"""
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            row.append(get_tsubame_image(scene_id, zoom, topleft_x + x, topleft_y + y))
        img.append(np.hstack(row))
    return  np.vstack(img)

def main():
    # 新国立競技場
    zoom=18
    lat, log = 35.6796057, 139.711851
    xtile, ytile = get_tile_num(lat, log, zoom)
    
    tsubame_image = get_tsubame_image(tsubame_scenes[0]["entityId"], zoom, xtile, ytile)

    img_4x4 = get_slat_series_image(tsubame_scenes[0]["entityId"], zoom, xtile, ytile, 4, 4)
    io.imshow(img_4x4)

if __name__ == '__main__':
    main()

秩父雲海発生予想エリア by ASTER GDEM 2

目的

秩父の雲海は標高の低いところに滞留すると仮定し、ASTER GDEM2 の標高モデルより、秩父雲海の発生エリアを推定する。

Tellus

https://www.tellusxdp.com/ja/

Result

秩父市標高データ　ASTER GDEM提供：METI and NASA
秩父市　部分的雲海の発生予想エリア
ASTER GDEM提供：METI and NASA
秩父市　全面的雲海の発生予想エリア
ASTER GDEM提供：METI and NASA

Code

import requests, json, math
from skimage import io
from io import BytesIO
import numpy as np

TOKEN = 'Enter your code'

# 秩父市
# 35.9827619,138.8036913
# 秩父市　ミューズパーク
zoom = 12
lat, lon = 35.9951626,139.0531834
y_slice = 160

def get_tile_num(lat_deg, lon_deg, zoom):
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    print("Zoom =", zoom, ", xtile =", xtile, ", ytile =", ytile)
    return (xtile, ytile)

def get_astergdem2_dsm(zoom, xtile, ytile):
    """
    """
    url = " https://gisapi.tellusxdp.com/astergdem2/dsm/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

def get_height_data(astergdem2_dsm):
    R = 255 * 255 * (astergdem2_dsm[:,:,0].astype(np.uint32))
    G = 255 * (astergdem2_dsm[:,:,1].astype(np.uint32))    
    B = astergdem2_dsm[:,:,2].astype(np.uint32)
    height =  R + G + B
    return height.astype(np.uint32)

def get_max_height(imgin):
    max_height = np.max(imgin)
    print("標高:", max_height)
    return max_height

def get_height_img(imgin, sea_level):
    gray =  imgin * 255.0 / get_max_height(imgin)
    _img = np.zeros([256,256,3],np.uint8)
    _img[:,:,0] = gray[:,:]
    _img[:,:,1] = gray[:,:]
    _img[:,:,2] = gray[:,:]

    # 標高 0m をブルーで塗りつぶす
    _img[:,:,2][_img[:,:,2] <= sea_level] = 255
    
    return _img.astype(np.uint8)

def get_height(zoom, lat,lon):
    (xtile, ytile) = get_tile_num(lat,lon, zoom)

    _astergdem2_dsm = get_astergdem2_dsm(zoom, xtile, ytile)

    _height = get_height_data(_img)
    
    return _height_data


def show_slice(height, y_slice):
    '''
    断面図を表示
    '''
    import matplotlib.pyplot as plt

    plt.figure(figsize=(10, 2.0))

    x = range(0,255)
    # 断面の Y座標

    plt.plot(x, height[y_slice, x].astype(np.float), label='mix', color='red')

    plt.grid(True)

    plt.show()
    
    
def main(zoom, lat, lon, y_slice, sea_level):
    (xtile, ytile) = get_tile_num(lat,lon, zoom)
    astergdem2_dsm = get_astergdem2_dsm(zoom, xtile, ytile)
    height_data = get_height_data(astergdem2_dsm)
    height_img = get_height_img(height_data, sea_level)
    io.imshow(height_img)
    show_slice(height_data, y_slice)

if __name__ == '__main__':
    main(zoom, lat, lon, y_slice, sea_level)

# 秩父市
# 35.9827619,138.8036913
# 秩父市　ミューズパーク
zoom = 12
lat, lon = 35.9951626,139.0531834
sea_level = 80

main(zoom, lat, lon, y_slice, sea_level)

ASTER GDEM 2

Tellus

https://www.tellusxdp.com/ja/

Result

ASTER GDEM提供：METI and NASA
ASTER GDEM提供：METI and NASA
ASTER GDEM提供：METI and NASA

Code

import requests, json, math
from skimage import io
from io import BytesIO
import numpy as np

TOKEN = 'Enter your token'

# 阿蘇山　32.8985056,131.0787207
# zoom = 9
# lat,lon = 32.8985056,131.0787207
# y_slice = 110

# 霧島山 31.9341693,130.8527772
# zoom = 9
# lat,lon = 31.9341693,130.8527772
# y_slice = 10

# 桜島 31.5833323,130.6412453
# zoom = 9
# lat,lon = 31.5833323,130.6412453
# y_slice = 155

# # イエローストーン 44.4122632,-110.7319387
# zoom = 8
# lat,lon = 44.4122632,-110.7319387
# y_slice = 110

# 富士山 35.3586732,138.719563
zoom = 9
lat,lon = 35.3586732,138.719563
y_slice = 44

def get_tile_num(lat_deg, lon_deg, zoom):
    # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    print("Zoom =", zoom, ", xtile =", xtile, ", ytile =", ytile)
    return (xtile, ytile)

def get_astergdem2_dsm(zoom, xtile, ytile):
    """
    """
    url = " https://gisapi.tellusxdp.com/astergdem2/dsm/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

def get_height_data(astergdem2_dsm):
    R = 255 * 255 * (astergdem2_dsm[:,:,0].astype(np.uint32))
    G = 255 * (astergdem2_dsm[:,:,1].astype(np.uint32))
    B = (astergdem2_dsm[:,:,2].astype(np.uint32))
    height =  R + G + B
    return height.astype(np.uint32)

def get_max_height(imgin):
    max_height = np.max(imgin)
    print("標高:", max_height)
    return max_height

def get_height_img(imgin):
    gray =  imgin * 255.0 / get_max_height(imgin)
    _img = np.zeros([256,256,3],np.uint8)
    _img[:,:,0] = gray[:,:]
    _img[:,:,1] = gray[:,:]
    _img[:,:,2] = gray[:,:]

    # 標高 0m をブルーで塗りつぶす
    _img[:,:,2][_img[:,:,2] == 0] = 255
    
    return _img.astype(np.uint8)

def get_height(zoom, lat,lon):
    (xtile, ytile) = get_tile_num(lat,lon, zoom)

    _astergdem2_dsm = get_astergdem2_dsm(zoom, xtile, ytile)

    _height = get_height_data(_img)
    
    return _height_data


def show_slice(height, y_slice):
    '''
    断面図を表示
    '''
    import matplotlib.pyplot as plt

    plt.figure(figsize=(10, 2.0))

    x = range(0,255)
    # 断面の Y座標

    plt.plot(x, height[y_slice, x].astype(np.float), label='mix', color='red')

    plt.grid(True)

    plt.show()
    
if __name__ == '__main__':
    main(zoom, lat, lon, y_slice)
    
def main(zoom, lat, lon, y_slice):
    (xtile, ytile) = get_tile_num(lat,lon, zoom)
    astergdem2_dsm = get_astergdem2_dsm(zoom, xtile, ytile)
    height_data = get_height_data(astergdem2_dsm)
    height_img = get_height_img(height_data)
    io.imshow(height_img)
    show_slice(height_data, y_slice)