2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

衛星データから東京23区の緑地率を算出し、それぞれの区の平均年所得との相関があるのか、調査をしてみました。

というのも、「東京23区における緑化の現状に関する研究」という研究から、東京23区それぞれの緑地域に住んでいる人の年収と緑地が近くにない人の年収のデータを見ると相関があったとの言葉があったのを見つけたから。それならば23区全体に範囲を拡げてもその傾向はそのままなのか？と気になり、検証を行った次第です。

本記事では、Tellusを活用して、東京23区を分割したうえで衛星データの画像を処理する方法を紹介しています。ぜひ読者のみなさんも衛星データから分かる東京23区それぞれの違いを様々な角度から比較、分析して遊んでみてください。

※本記事は宙畑メンバーが気になったヒト・モノ・コトを衛星画像から探す不定期連載「宇宙データ使ってみた-Space Data Utilization-」の第17弾です。まだまだ修行中の身のため至らない点があるかと思いますがご容赦・アドバイスいただけますと幸いです！

(1)論文の内容を確認する

まず、「東京23区における緑化の現状に関する研究」について、その内容を簡単にまとめると以下のような内容が記してありました。
「墨田区を除き、東京23区それぞれの緑化率の高い地域と低い地域において、その緑化率と、野村総合研究所のエリア別の推計結果「世帯当たり年間総所得」を参考にした（緑被率の高い地域，低い地域）における年間総所得を並べた時に相関が見られる」

そこで、「空港周辺の緑地率とバードストライク率の相関を検証してみた」にて緑地率を算出した方法を応用して、実際にそれを確かめてみようという訳です。

(2)東京23区の平均年所得ランキングを見る

まず、東京23区の平均年収を確認しましょう。

「年収ガイド」によると、平成2018年の東京23区別平均所得は以下のようになっておりました。
※各種の社会保険や控除などがあるため実際の総年収としては、この「所得」数値よりも多くなるとのこと

AVNIR2の植生データ(緑色が植物のある場所)を表示したTellus上で東京の行政区分のgeojsonを読み込んだものが以下。

目視で見比べてみると、渋谷区や港区といった大きな公園がある区はたしかに年収が高いように見えます。

では、実際にはどうなのでしょうか。実際に衛星データから緑地率を算出して確かめてみましょう。

(3)データの説明

今回も、「空港周辺の緑地率とバードストライク率の相関を検証してみた」と同じく、Tellusに搭載されているAVNIR-2という光学センサのデータを活用して、東京23区の緑地率を調べました。

衛星データでなぜ緑地率を見れるのか、詳しくは「光の波長って何？ なぜ人工衛星は人間の目に見えないものが見えるのか」をご覧ください。

また、今回は東京23区の区域に絞って緑地率を算出するということが新しいチャレンジ。そのコードの解説と合わせて、以下ご覧ください。

(4)コード解説&取得

それでは実際に、東京23区それぞれの緑地率を確認していきましょう。

まずは必要なライブラリをインストールしておきます。

import requests, math
import numpy as np
from pathlib import Path
from skimage import io
from skimage.draw import polygon
from io import BytesIO
%matplotlib inline

TOKEN = '自分のトークンを記載'

次に、手元にある東京23区のgeojsonに対応する衛星データの情報とタイル情報を確認し、画像で吐き出すためのコードを書きます。

tokyo_wards = [{
    'name': 'chiyoda',
    'z': 14,
    'x': 14551,
    'y': 6450,
    'size_x': 3,
    'size_y': 3,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13101'
}, {
    'name': 'chuoku',
    'z': 14,
    'x': 14552,
    'y': 6450,
    'size_x': 2,
    'size_y': 3,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13102'
}, {
    'name': 'minatoku',
    'z': 14,
    'x': 14550,
    'y': 6451,
    'size_x': 2,
    'size_y': 3,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13103'
}, {
    'name': 'shinagawaku',
    'z': 14,
    'x': 14549,
    'y': 6451,
    'size_x': 4,
    'size_y': 4,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13109'
}, {
    'name': 'ohtaku',
    'z': 14,
    'x': 14547,
    'y': 6455,
    'size_x': 7,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13111'
}, {
    'name': 'setagayaku',
    'z': 14,
    'x': 14544,
    'y': 6451,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152880',
    'rspId': 'D072P0',
    'code': '13112'
}, {
    'name': 'suginamiku',
    'z': 14,
    'x': 14544,
    'y': 6448,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13115'
}, {
    'name': 'nakanoku',
    'z': 14,
    'x': 14546,
    'y': 6448,
    'size_x': 4,
    'size_y': 4,
    'productId': 'ALAV2A028152880',
    'rspId': 'D072P0',
    'code': '13114'
}, {
    'name': 'shinjukuku',
    'z': 14,
    'x': 14548,
    'y': 6448,
    'size_x': 4,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13104'
}, {
    'name': 'toshimaku',
    'z': 14,
    'x': 14548,
    'y': 6448,
    'size_x': 4,
    'size_y': 2,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13116'
}, {
    'name': 'nerimaku',
    'z': 14,
    'x': 14543,
    'y': 6446,
    'size_x': 6,
    'size_y': 4,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13120'
}, {
    'name': 'itabashiku',
    'z': 14,
    'x': 14546,
    'y': 6444,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13119'
}, {
    'name': 'kitaku',
    'z': 14,
    'x': 14549,
    'y': 6445,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13117'
}, {
    'name': 'bunkyoku',
    'z': 14,
    'x': 14550,
    'y': 6448,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13105'
}, {
    'name': 'daitoku',
    'z': 14,
    'x': 14552,
    'y': 6448,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13106'
}, {
    'name': 'sumidaku',
    'z': 14,
    'x': 14553,
    'y': 6448,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13107'
}, {
    'name': 'kotoku',
    'z': 14,
    'x': 14553,
    'y': 6450,
    'size_x': 5,
    'size_y': 7,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13108'
}, {
    'name': 'edogawaku',
    'z': 14,
    'x': 14555,
    'y': 6447,
    'size_x': 5,
    'size_y': 9,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13123'
}, {
    'name': 'katsushikaku',
    'z': 14,
    'x': 14555,
    'y': 6445,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13122'
}, {
    'name': 'adachiku',
    'z': 14,
    'x': 14551,
    'y': 6444,
    'size_x': 7,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13121'
}, {
    'name': 'arakawaku',
    'z': 14,
    'x': 14552,
    'y': 6447,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13118'
}, {
    'name': 'meguroku',
    'z': 14,
    'x': 14546,
    'y': 6452,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13110'
}, {
    'name': 'shibuyaku',
    'z': 14,
    'x': 14548,
    'y': 6451,
    'size_x': 5,
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13113'
}]


def get_combined_image(get_image_func, z, topleft_x, topleft_y, size_x=1, size_y=1, option={}, params={}):
    rows = []
    blank = np.zeros((256, 256, 4), dtype=np.uint8)

    for y in range(size_y):
        row = []
        for x in range(size_x):
            try:
                img = get_image_func(z, topleft_x + x, topleft_y + y, option, params)
            except Exception as e:
                img = blank
                
            row.append(img)
        rows.append(np.hstack(row))
    return  np.vstack(rows)


def get_AVNIR2_scene_tile(zoom, xtile, ytile, option, params={}):
    url = 'https://gisapi.tellusxdp.com/av2ori/{}/{}/{}/{}/{}.png'.format(option['rspId'], option['productId'], zoom, xtile, ytile)
    headers = {'Authorization': 'Bearer ' + TOKEN}

    r = requests.get(url, params=params, headers=headers)

    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return io.imread(BytesIO(r.content))

tokyo_imgs = {}
for word in tokyo_wards:
    # タイル取得
    img = get_combined_image(get_AVNIR2_scene_tile, word['z'], word['x'], word['y'],
                         word['size_x'], word['size_y'], 
                         option={'productId':word['productId'],'rspId':word['rspId']}, params={'preset': 'ndvi'})  
    tokyo_imgs[word['name']] = img

今回は東京23区がすっぽりと収まる衛星データが見つかったため、そのシーンは固定で使用し、ズームレベルを14に統一し、タイルを連結させて23区それぞれのデータに対応するタイル情報をまとめました。
※シーンを統一しなければ、シーンによっては色味も変わってしまい、緑地率の計算もくるってしまうため、要注意です
※ズームレベルを合わせない場合も、1枚の画像の解像度が変化することで色味が変わってしまうため、緑地率の計算もくるってしまうため、要注意です

たとえば、新宿区の場合は、左上から右に4枚、下に5枚のタイルを連結させると新宿区がすっぽり収まる連結タイルが作成できるため、以下のようになります。

各区のタイル連結例として新宿区を抜粋

{
    'name': 'shinjukuku',
    'z': 14,
    'x': 14548, # 左上のx座標
    'y': 6448,
    'size_x': 4, # 右に●枚
    'size_y': 5,
    'productId': 'ALAV2A028152870',
    'rspId': 'D072P0',
    'code': '13104'
}

上記のようにすべての区のタイル情報を記載終えたら次のコードを書きましょう。次は、geojsonファイルの読み込みと、それぞれの区のクリッピングを行います。

def find_features(filename, code):
    with open(filename) as f:
        df = json.load(f)
        return [features for features in df['features'] if features['properties']['N03_007']==code]
    return []

def get_tile_bbox(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """
    タイル座標からバウンダリボックスを取得
    https://tools.ietf.org/html/rfc7946#section-5
    """
    def num2deg(xtile, ytile, zoom):
        # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        return (lon_deg, lat_deg)
    
    right_top = num2deg(topleft_x + size_x , topleft_y, zoom)
    left_bottom = num2deg(topleft_x, topleft_y + size_y, zoom)
    return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])

def get_polygon_image(points, bbox, imgsize):
    def world2pixel(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(world2pixel(bbox, (imgsize['width'], imgsize['height']), p[0], p[1]))
    pixels = np.array(pixels)

    poly = np.ones((imgsize['height'], imgsize['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, pxsize):
    polygon_imgs = []
    for feature in features:
        points = feature['geometry']['coordinates'][0]
        polygon_imgs.append(get_polygon_image(points, bbox, pxsize))

    mask = np.where((polygon_imgs[0] > 0), 0, 1)
    for i in range(1,len(polygon_imgs)):
        mask += np.where((polygon_imgs[i] > 0), 0, 1)

    return np.where((mask > 0), 0, 1).astype(np.uint8)


def clip(img, mask):
    base = img.transpose(2, 0, 1)[:3,:]
    cliped = np.choose(mask, (base, 0)).astype(np.uint8)
    return cliped.transpose(1, 2, 0)


for word in tokyo_wards:
    # geojson読み込む
    features = find_features('N03-19_13_190101.geojson', word['code'])
    
    # クリッピング
    pxsize = {'width': 256*word['size_x'], 'height': 256*word['size_y']}
    bbox = get_tile_bbox(word['z'], word['x'], word['y'], word['size_x'], word['size_y'])
    mask = get_mask_image(features, bbox, pxsize)
    tokyo_imgs[word['name']] = clip(tokyo_imgs[word['name']], mask)

io.imshow(tokyo_imgs['kotoku'])

ここでは、まず指定した連結タイルのピクセル情報と4隅の緯度経度情報を取得して、各区がすっぽりと入るように連結したタイルのサイズを把握します。その後、東京23区それぞれのポリゴン情報(features)からその特徴を構成する座標を把握して、連結タイルの中にその情報を読み込み、それぞれの区のデータをクリップする記述をしています。「features」と複数形になっているのは、飛び地でもしっかりと同じ区として取得できるようになっているためです。

また、実際に飛び地のある江東区を確認する記述も記載しており、本当にそのようなコードとなっているかを確かめた結果出てきた画像がこちら。

上記のように、飛び地があってもきれいに切り取れていますね。

それではいよいよ緑地率の算出を行います。まずはコードから。

import cv2

def count_not_black(img):
    return len(np.nonzero(img.ravel())[0])

def filter_by_color(img, color_min, color_max):
    # OpenCVで扱うために色順をRGB(RGBA)からBGRに入れ替える
    cv_bgr = cv2.cvtColor(img, cv2.COLOR_RGBA2BGR) 
    # BGRからHSVに変換
    cv_hsv = cv2.cvtColor(cv_bgr, cv2.COLOR_BGR2HSV)
    
    # 指定した範囲の色以外を黒に変換する
    mask = cv2.inRange(cv_hsv, color_min, color_max)
    masked_img = cv2.bitwise_and(cv_bgr, cv_bgr, mask=mask)

    # 得られた画像の色順をBGRからRGBに戻す
    return cv2.cvtColor(masked_img, cv2.COLOR_BGR2RGB)


def calc_vegetation_rate(img):
    green_min = np.array([25,0,0]) 
    green_max = np.array([70,255,255])
    filtered = filter_by_color(img, green_min, green_max)
    return count_not_black(filtered)/count_not_black(img)


for k, v in tokyo_imgs.items():
    # 植被率
    print(k + ' : ' + str(calc_vegetation_rate(v)))

ここでは、「空港周辺の緑地率とバードストライク率の相関を検証してみた」と同様、画像を変換したうえで、それぞれの区で切り取ったすべての画像から、指定の閾値で緑地と判断した部分の割合を返すコードとなっており、東京23区の緑地率を取得することができます。

これで東京23区別の緑地率と平均年収が出そろいました。いよいよ相関を検証してみます。

(5)結果

さっそく、相関を検証してみた結果が以下のようになります。

相関係数：0.658
P値：0.001

やや相関はあるようですが、ばらつきはあるようです。ここからさらに、下町や練馬区など、いわゆる都心と呼ばれる場所からは少し外れる区を外した都心16区と呼ばれる区のみに絞るとどうかなと思い、再度相関を検証してみた結果がこちら。

相関係数：0.716
P値：0.002

相関係数の値がやや大きくなりました。

他にも、緑地率が高いということは公園が多い？つまり、子供を育てやすい環境にあるのでは？ということで、東京23区それぞれの出生率と緑地率の相関を検証してみた結果がこちらです。

相関係数：0.117
P値：0.596

結果を見ると、ほとんど関係なさそうですね。

では、緑地が多いところでは犬の散歩がしやすいということで、世帯に対する犬の登録頭数の割合を犬の所有率として出したときに相関があるのでは？と調べてみたところ、以下のようになりました。

相関係数：0.267
P値：0.317

こちらも相関があります！とは断言できない結果となりました。公園付近など、より地域を絞るとまた違った結果がなるのかもしれません。

(6)まとめ

東京23区の緑地率と年収に相関はあるか？と思って検証してみた本記事でしたが、冒頭の研究結果のような目を見張るような結果を得ることはできませんでした。参考にした研究では、東京23区の中からさらに緑地率の高い地域と低い地域を絞っての比較をしていたため、年収に限らず今回例として挙げた項目との相関を見ても、区という行政区分ではいささか広過ぎるのかもしれません。

今回は東京23区をクリッピングして、緑地率を算出する方法をご紹介しましたが、応用すれば水の割合やその他気になるデータも同様の手順で取得できます。

例えば、今の時期はインフルエンザの感染のしやすさや、家賃との相関なども、衛星データから得られる様々な情報を組み合わせて調べられると、興味深い検証結果を得られるかもしれません。

皆さんも東京23区やもっと細かい行政区分ごとなどで違いを比較したいといった何か気になることがあれば、Tellusを活用してみませんか？

Tellusの登録はこちらから