Tellus画像処理の練習メモ：緑地を切り出してGeoJSONを生成してみる

Tellus-APiの方はだいたい理解できたので、次はJupyterNoteでの画像処理に挑戦してみました。
Pythonでの画像処理は初めてなので、pashango2さんのエントリ「PIL/Pillow チートシート」を参考にしました。なお、緑地抽出に適した処理手順はわからないので、色々と施行錯誤しながら緑地抽出を試みてみました。

TellusOSを使って、衛星画像の画像処理を練習中です。光学画像から緑地を抽出し、GeoJSONを生成する処理を試してみました。https://t.co/QZr92WxBAQ#Tellus pic.twitter.com/3xqLna9jTW

— t_mat (@t_mat) April 1, 2019

##画像処理

####1.元画像の読み込み
ASNARO-01「東京北部2016/05/05」から、適度に緑地があるが場所（東京メトロ王子駅付近）をスクリーンショットし、その画像（/tellus_data/に保存されている）をJupyterNoteに読み込みました。

from PIL import Image, ImageFilter ,ImageOps,ImageMath,ImageChops
import numpy as np
import matplotlib.pyplot as plt
im = Image.open('../../tellus_data/20190331-222445.tiff','r')
plt.imshow(np.array(im))
im.save('../data/20190331-222445_0.png')

####2.画素の強調
元画像がくすんでいるので、同じ画像を4回ソフトライトして画素を強調してみました。

def _blend_f(img1, img2, func):
    blend_eval = "convert(func(float(a), float(b)), 'L')"
    bands = [
        ImageMath.eval(
            blend_eval,
            a=a,
            b=b,
            func=func
        )
        for a, b in zip(img1.split(), img2.split())
    ]
    return Image.merge(img1.mode, bands)

def _soft_light(a, b):
    _cl = (a / 255) ** ((255 - b) / 128) * 255
    _ch = (a / 255) ** (128 / b) * 255
    return _cl * (b < 128) + _ch * (b >= 128)
im2=_blend_f(im, im, _soft_light)
plt.imshow(np.array(im2))
im3=_blend_f(im2, im2, _soft_light)
plt.imshow(np.array(im3))

####3.ポスタライズ
ポスタライズで色を単純化しました。

im4=im3.convert("RGB")
im5=ImageOps.posterize(im4, 4)
plt.imshow(np.array(im5))
im5.save('../data/20190331-222445_2.png')

####4.ガンマ補正
ガンマ補正で照度、彩度を調整しました。

def gamma_table(gamma_r, gamma_g, gamma_b, gain_r=1.0, gain_g=1.0, gain_b=1.0):
    r_tbl = [min(255, int((x / 255.) ** (1. / gamma_r) * gain_r * 255.)) for x in range(256)]
    g_tbl = [min(255, int((x / 255.) ** (1. / gamma_g) * gain_g * 255.)) for x in range(256)]
    b_tbl = [min(255, int((x / 255.) ** (1. / gamma_b) * gain_b * 255.)) for x in range(256)]
    return r_tbl + g_tbl + b_tbl

im6=im5.point(lambda x: x * 2.0)
plt.imshow(np.array(im6))
im6.save('../data/20190331-222445_3.png')

####5.画素の強調
もう一度ソフトライトを行いました。

im7=_blend_f(im6, im6, _soft_light)
plt.imshow(np.array(im7))
im7.save('../data/20190331-222445_4.png')

####6.色抽出その1
OpenCVをロードし、緑系の色のピクセルを抽出しました。

!pip install opencv-python
import cv2
image_file = '../data/20190331-222445_4.png'
im8 = cv2.imread(image_file)
hsv = cv2.cvtColor(im8, cv2.COLOR_BGR2HSV)
lower = np.array([32, 64, 32])
upper = np.array([180, 255, 120])
img_mask = cv2.inRange(hsv, lower, upper)
im9 = cv2.bitwise_and(im8, im8, mask=img_mask)
plt.imshow(np.array(im9))
cv2.imwrite('../data/20190331-222445_5.png', im9)

####7.色抽出その2
上記では黄緑系統の色が抽出されなかったので、再度、色抽出を行いました。

lower = np.array([32, 64, 32])
upper = np.array([200, 255, 120])
img_mask = cv2.inRange(im8, lower, upper)
im10 = cv2.bitwise_and(im8, im8, mask=img_mask)
plt.imshow(np.array(im10))
cv2.imwrite('../data/20190331-222445_6.png', im10)

####8.色抽出結果の合成
6、7で抽出した画像を合成しました。

im11 = Image.open('../data/20190331-222445_5.png','r')
im12 = Image.open('../data/20190331-222445_6.png','r')
im13=ImageChops.lighter(im11, im12)
plt.imshow(np.array(im13))
im13.save('../data/20190331-222445_7.png')

####9.ピクセルの再抽出
画像のピクセルを検査し、RGBのうちG値が最も高いピクセルを抽出しました。

size=im13.size
im14 = Image.new('RGB',size)
for x in range(size[0]):
    for y in range(size[1]):
        r,g,b = im13.getpixel((x,y))
        if( g > r and g > b):
            im14.putpixel((x,y),(r,g,b))
        else:
            im14.putpixel((x,y),(0,0,0))
plt.imshow(np.array(im14))
im14.save('../data/20190331-222445_8.png')

####10.収縮・膨張
収縮・膨張により、ピクセル上のごみを処理し、その後、ガンマ補正を行いました。

im15=im14.filter(ImageFilter.MinFilter()) 
im16=im15.filter(ImageFilter.MaxFilter()) 
plt.imshow(np.array(im16))
im16.save('../data/20190331-222445_9.png')
im17=im16.point(lambda x: x * 2.0)
plt.imshow(np.array(im17))
im17.save('../data/20190331-222445_10.png')

####11.減色処理
色を黒を含め、4色に減色処理しました。

def sub_color(src, K):
    Z = src.reshape((-1,3))
    Z = np.float32(Z)
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)
    ret, label, center = cv2.kmeans(Z, K, None, criteria, 10, cv2.KMEANS_RANDOM_CENTERS)
    center = np.uint8(center)
    res = center[label.flatten()]
    return res.reshape((src.shape))
im19 = cv2.imread('../data/20190331-222445_11.png') 
im20 = sub_color(im19, 4)
plt.imshow(np.array(im20))
cv2.imwrite('../data/20190331-222445_12.png', im20)

####12.色の分離
11で緑3色に減色したので、各色のピクセルを抽出して分離し、2値画像に変換しました。
なお、この処理は使い慣れたJavaで実施しました。

####13.ベクター化
上記の2値画像をpotraceでSVGに変換しました。

####14.SVG->GeoJSON生成
svg2geojsonにより、SVGをGeoJSONに変換しました。
Tellusに生成したGeoJSONを読み込ませた画面を示しますが、そこそこうまく抽出できています。ただ、感覚的ですが、5％ほど、緑地の抽出漏れがある感じです。

##雑感

• 光学画像からでもそこそこの精度で緑地が抽出でき、GeoJSONを生成できることがわかりました。また、林相と画像の色相に関係性があるようなら、林相区分等の自動処理も可能かもと思いました。
• 自動処理で緑地を抽出するとなると、処理手順や抽出の閾値を色々と検討・工夫する必要があると感じました。
• ベクター化やGeoJSON生成には外部ライブラリを使用しましたが、独自に処理を実装するか機能拡張した方がスムーズかもと感じました。