任意地点の標高データを取得・解析する

Tag

世田谷区の標高のばらつきを算出する

グラフを作る以外にも標高データを使って学んでみましょう。

各区の形状は国土数値情報ダウンロードサービスの行政区域データを利用します。行政区域データのGeoJSONのダウンロード方法はこちらの記事を参考にしてください。

本記事では「work」ディレクトリの下に「data」ディレクトリを作り、そこにGeoJSONをアップロードしています。

devAster_proj_20200220_1.png

このGeoJSONを読み込み、中から自治体コードが13112の要素を抜き出します。

import os
from pathlib import Path
basepath = Path(os.path.expanduser('~'))
filepath = basepath / 'work/data/N03-19_13_190101.geojson'
# ダウンロードした行政区域データを呼び出す
with open(filepath) as f:
    df = json.load(f)
# 世田谷区(市区町村コード13112)のfeatureを抜き出す
tokyo_13112_features = [f for f in df['features'] if f['properties']['N03_007']=='13112']
print(len(tokyo_13112_features))
print(tokyo_13112_features[0])

devAster_proj_20200220_2.png 続いて、世田谷区全域が収まるサイズで標高タイルを取得します。ズーム12で縦3枚、横2枚の計6枚をつなぎ合わせます。

import numpy as np
def get_connected_astergdem2_dsm(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """
    ASTER GDEM2 の標高タイルを size_x * size_y 分つなげた1枚の画像を作る
    """ 
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            row.append(get_astergdem2_dsm(zoom, topleft_x + x, topleft_y + y))
        img.append(np.hstack(row))
    return  np.vstack(img)
tokyo_13112_img = get_connected_astergdem2_dsm(12, 3636, 1612, 2, 3)
io.imshow(tokyo_13112_img)

devAster_proj_20200220_4.png

標高が全体的に低いことから、全体的に黒っぽい画像となっています
ASTER GDEM提供:METI and NASA

次に世田谷区のポリゴンから画像作成し、それを重ね合わせることで世田谷区の形状に切り抜かれた標高タイルを作成します。

from skimage.draw import polygon
def get_polygon_image(points, bbox, img_size):
    """
    bboxの領域内にpointsの形状を描画した画像を作成する
    """ 
    def lonlat2pixel(bbox, size, lon, lat):
        """
        経度緯度をピクセル座標に変換
        """
        dist = ((bbox[2] - bbox[0])/size[0], (bbox[3] - bbox[1])/size[1])
        pixel = (int((lon - bbox[0]) / dist[0]), int((bbox[3] - lat) / dist[1]))
        return pixel
    pixels = []
    for p in points:
        pixels.append(lonlat2pixel(bbox, (img_size['width'], img_size['height']), p[0], p[1]))
    pixels = np.array(pixels)
    poly = np.ones((img_size['height'], img_size['width']), dtype=np.uint8)
    rr, cc = polygon(pixels[:, 1], pixels[:, 0], poly.shape)
    poly[rr, cc] = 0
    return poly
def get_mask_image(features, bbox, px_size):
    """
    図形からマスク画像を作成する
    Parameters
    ----------
    features : array of GeoJSON feature object 
        マスクするfeature objectの配列 
    bbox : array
        描画領域のバウンディングボックス 
    px_size : array 
        描画領域のピクセルサイズ
    Returns
    -------
    img_array : ndarray
        マスク画像
    """
    ims = []
    for i in range(len(features)):
        points = features[i]['geometry']['coordinates'][0]
        ims.append(get_polygon_image(points, bbox, px_size))
    mask = np.where((ims[0] > 0), 0, 1)
    for i in range(1,len(ims)):
        mask += np.where((ims[i] > 0), 0, 1)
    return np.where((mask > 0), 0, 1).astype(np.uint8)
bbox = get_tile_bbox(12, 3636, 1612, 2, 3)
px_size = {'width': tokyo_13112_img.shape[1], 'height': tokyo_13112_img.shape[0]}
tokyo_13112_mask = get_mask_image(tokyo_13112_features, bbox, px_size)
io.imshow(tokyo_13112_mask)

devAster_proj_20200220_3.png

ASTER GDEM提供:METI and NASA
def clip_image(img, mask):
    """
    図形をマスクする
    Parameters
    ----------
    img : ndarray
        加工される画像 
    mask : ndarray
        マスク画像 
    Returns
    -------
    img_array : ndarray
        マスクされた画像
    """
    base = img.transpose(2, 0, 1)[:3,:]
    cliped = np.choose(mask, (base, 0)).astype(np.uint8)
    return cliped.transpose(1, 2, 0)
tokyo_13112_cliped = clip_image(tokyo_13112_img, tokyo_13112_mask)

区の形状に切り抜かれた標高タイルを利用してヒストグラムを作成します。

def get_height_array(tile, u):
    """
    標高タイルから高さの配列を作る
    ただしエラーや0メートル地点は配列から除く
    """
    heights = []
    heights_err = []
    for i in range(tile.shape[0]):
        for j in range(tile.shape[1]):
            R = tile[i, j, 0]
            G = tile[i, j, 1]
            B = tile[i, j, 2]
            try:
                heights.append(calc_height_chiriin_style(R, G, B, u))
            except ValueError as e:
                heights.append(0)
                heights_err.append((i,j))
    heights = np.array(heights)
    return (heights[heights > 0], heights_err)  
def plot_height_hist(tile, u=1):
    """
    標高タイルからヒストグラムを作る
    Parameters
    ----------
    tile : ndarray 
        標高タイル
    """
    (heights, heights_err) = get_height_array(tile, u)
    fig, ax= plt.subplots(figsize=(8, 4)) 
    plt.hist(heights, bins=30, rwidth=0.95)
    ax.text(0.01, 0.95, 'Errors: {}'.format(len(heights_err)), transform=ax.transAxes)
    plt.show()
plot_height_hist(tokyo_13112_cliped)

devAster_proj_20200220_5.png

ASTER GDEM提供:METI and NASA

ヒストグラムを見ますと、40から50mのあたりに値が集中していることがわかります。分布はやや右に長い裾を引いており、度数自体は左に偏りがあり、いわゆる正規分布からは外れた分布だとわかります。

最後に分散や平均も求めてみましょう。

def calc_height_statistics(tile, u=1):
    """
    標高タイルの高さの統計情報を求める
    Parameters
    ----------
    tile : ndarray 
        標高タイル 
    Returns
    -------
    result_object : object
        統計情報
    """
    (heights, heights_err) = get_height_array(tile, u)
    return {'std': np.std(heights), 'mean': np.mean(heights), 'max': np.max(heights), 'min': np.min(heights)}
print(calc_height_statistics(tokyo_13112_cliped))

devAster_proj_20200220_6.png

他の区についても同じ手法で求めることができます。

以上が、TellusのJupyterLabで富士山の標高グラフを作成する方法をご紹介しました。ASTER GDEM 2.0のデータは、国土地理院が提供する標高データに精度で劣りますが、データが欠落している範囲が少なく、国外の地点もカバーできている点で優れています。ぜひ活用してください。

今回使用したスクリプトです。

import os
from pathlib import Path
basepath = Path(os.path.expanduser('~'))
filepath = basepath / 'work/data/N03-19_13_190101.geojson'
# ダウンロードした行政区域データを呼び出す
with open(filepath) as f:
    df = json.load(f)
# 世田谷区(市区町村コード13112)のfeatureを抜き出す
tokyo_13112_features = [f for f in df['features'] if f['properties']['N03_007']=='13112']
print(len(tokyo_13112_features))
print(tokyo_13112_features[0])
import numpy as np
def get_connected_astergdem2_dsm(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """
    ASTER GDEM2 の標高タイルを size_x * size_y 分つなげた1枚の画像を作る
    """
    img = []
    for y in range(size_y):
        row = []
        for x in range(size_x):
            row.append(get_astergdem2_dsm(zoom, topleft_x + x, topleft_y + y))
        img.append(np.hstack(row))
    return  np.vstack(img)
tokyo_13112_img = get_connected_astergdem2_dsm(12, 3636, 1612, 2, 3)
io.imshow(tokyo_13112_img)
from skimage.draw import polygon
def get_polygon_image(points, bbox, img_size):
    """
    bboxの領域内にpointsの形状を描画した画像を作成する
    """
    def lonlat2pixel(bbox, size, lon, lat):
        """
        経度緯度をピクセル座標に変換
        """
        dist = ((bbox[2] - bbox[0])/size[0], (bbox[3] - bbox[1])/size[1])
        pixel = (int((lon - bbox[0]) / dist[0]), int((bbox[3] - lat) / dist[1]))
        return pixel
    pixels = []
    for p in points:
        pixels.append(lonlat2pixel(bbox, (img_size['width'], img_size['height']), p[0], p[1]))
    pixels = np.array(pixels)
    poly = np.ones((img_size['height'], img_size['width']), dtype=np.uint8)
    rr, cc = polygon(pixels[:, 1], pixels[:, 0], poly.shape)
    poly[rr, cc] = 0
    return poly
def get_mask_image(features, bbox, px_size):
    """
    図形からマスク画像を作成する
    Parameters
    ----------
    features : array of GeoJSON feature object
        マスクするfeature objectの配列
    bbox : array
        描画領域のバウンディングボックス
    px_size : array
        描画領域のピクセルサイズ
    Returns
    -------
    img_array : ndarray
        マスク画像
    """
    ims = []
    for i in range(len(features)):
        points = features[i]['geometry']['coordinates'][0]
        ims.append(get_polygon_image(points, bbox, px_size))
    mask = np.where((ims[0] > 0), 0, 1)
    for i in range(1,len(ims)):
        mask += np.where((ims[i] > 0), 0, 1)
    return np.where((mask > 0), 0, 1).astype(np.uint8)
bbox = get_tile_bbox(12, 3636, 1612, 2, 3)
px_size = {'width': tokyo_13112_img.shape[1], 'height': tokyo_13112_img.shape[0]}
tokyo_13112_mask = get_mask_image(tokyo_13112_features, bbox, px_size)
io.imshow(tokyo_13112_mask)
def clip_image(img, mask):
    """
    図形をマスクする
    Parameters
    ----------
    img : ndarray
        加工される画像
    mask : ndarray
        マスク画像
    Returns
    -------
    img_array : ndarray
        マスクされた画像
    """
    base = img.transpose(2, 0, 1)[:3,:]
    cliped = np.choose(mask, (base, 0)).astype(np.uint8)
    return cliped.transpose(1, 2, 0)
tokyo_13112_cliped = clip_image(tokyo_13112_img, tokyo_13112_mask)
def get_height_array(tile, u):
    """
    標高タイルから高さの配列を作る
    ただしエラーや0メートル地点は配列から除く
    """
    heights = []
    heights_err = []
    for i in range(tile.shape[0]):
        for j in range(tile.shape[1]):
            R = tile[i, j, 0]
            G = tile[i, j, 1]
            B = tile[i, j, 2]
            try:
                heights.append(calc_height_chiriin_style(R, G, B, u))
            except ValueError as e:
                heights.append(0)
                heights_err.append((i,j))
    heights = np.array(heights)
    return (heights[heights > 0], heights_err)  
def plot_height_hist(tile, u=1):
    """
    標高タイルから高さのヒストグラムを作る
    Parameters
    ----------
    tile : ndarray
        標高タイル
    """
    (heights, heights_err) = get_height_array(tile, u)
    fig, ax= plt.subplots(figsize=(8, 4))
    plt.hist(heights, bins=30, rwidth=0.95)
    ax.text(0.01, 0.95, 'Errors: {}'.format(len(heights_err)), transform=ax.transAxes)
    plt.show()
plot_height_hist(tokyo_13112_cliped)
def calc_height_statistics(tile, u=1):
    """
    標高タイルの高さの統計情報を求める
    Parameters
    ----------
    tile : ndarray
        標高タイル
    Returns
    -------
    result_object : object
        統計情報
    """
    (heights, heights_err) = get_height_array(tile, u)
    return {'std': np.std(heights), 'mean': np.mean(heights), 'max': np.max(heights), 'min': np.min(heights)}
print(calc_height_statistics(tokyo_13112_cliped))