Pythonで線形補間（欠けている標高データを補完してみる）

前回、みんなで使える無料の地図データOpenStreetMapと、フリーで利用可能な世界の陸域と海の地形データETOPO1グローバル地形データセットでコピーフリーの海底地形図付き日本地図を作ってみましたが

memomemokun.hateblo.jp

ETOPO1グローバル地形データセットの標高メッシュだと地図タイルのズームレベルで、ズームレベル６程度の解像度しか無いため、解像度の足りない部分は線形補間でそれらしい値を補いながら、ETOPO1の標高メッシュを拡張してみることに。

やっていることとしては、

netCDF4ライブラリーでETOPO1データを読み込み
日本周辺の地形図を描くのに必要な東経110度から東経170度、北緯16度から北緯56度の範囲でデータを切り出し（前回、データの切り出しまでfor文で２億回ほどグリグリ回してしまいましたが^^;; numpyの配列操作だとめちゃ早い）
切り出した標高データーを経度方向および緯度方向にそれぞれ１行１列だけメッシュの要素数を拡張しながら線形補間してみることに

まずは、以下ソースで実験。

#!/usr/bin/env python
# -*- coding: utf-8 -*- 

from math import pi
from math import tanh
from math import sin
from math import asin
from numpy import arctanh

from netCDF4 import Dataset
import numpy as np


# 経緯度のピクセル座標を求める
def fromLatLngToPoint(lat,lon, z):
    L = 85.05112878
    pointX = int( (2**(z+7)) * ((lon/180) + 1) )
    pointY = int( (2**(z+7)/pi) * (-1 * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180))) )
    return pointX, pointY


# 経緯度の世界座標を求める
def fromLatLngToWordLatLon(lat, lon, z):
    L = 85.05112878
    wordLon = (2**(z+7)) * ((lon/180) + 1)
    wordLat = (2**(z+7)/pi) * (-1. * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180)))
    return wordLat, wordLon
    

# ズームレベル
z = 7

#地形データを読み込む
etopo1 = Dataset('ETOPO1_Bed_g_gmt4.grd','r')

lons   = etopo1.variables['x'][:]
lats   = etopo1.variables['y'][:]
depths = etopo1.variables['z'][:]


print lats.shape
print lons.shape
print depths.shape

# 経緯度条件
condLons = (lons >= 110) * (lons <= 170)
condLats = (lats >= 16) * (lats <= 56)

# 経緯度条件を満たす配列インデックスを抽出
lonsIX = np.array(np.where(condLons)).flatten()
latsIX = np.array(np.where(condLats)).flatten()

# 経緯度条件を満たす経緯度と標高データを取り出す
lons = lons[lonsIX].flatten()
lats = lats[latsIX].flatten()
depths = depths[latsIX[0]:latsIX[-1]+1, lonsIX[0]:lonsIX[-1]+1]

print lats.shape
print lons.shape
print depths.shape

# 経緯度を世界座標に変換する
WordLats = np.array([fromLatLngToWordLatLon(lat,0, z)[0] for lat in lats])
WordLons = np.array([fromLatLngToWordLatLon(0,lon, z)[1] for lon in lons])

xdiff =  int(WordLons[-1]) - int(WordLons[0])
ydiff =  int(WordLats[0]) - int(WordLats[-1])


from scipy import interpolate
from matplotlib import pylab as plt

# 拡張したい経度方向、緯度方向の座標を生成
# 補完後の座標がピクセル座標ピッタリになるよう間隔は1づつの整数値になるようにする
pointX = np.linspace(int(WordLons[0]), int(WordLons[-1]), xdiff+1)
pointY = np.linspace(int(WordLats[0]), int(WordLats[-1]), ydiff+1)

# 標高を経度方向に線形補間しながら拡張する
depthsX = depths[0,:]
linearf = interpolate.interp1d(WordLons, depthsX, kind='linear', fill_value='extrapolate')
depthsExtendX = linearf(pointX)

plt.figure(figsize=(12, 9))
plt.plot(WordLons, depthsX, "o", label="original")
plt.plot(pointX, depthsExtendX, "r", label="linear")

plt.legend()
plt.savefig('1.jpg', dpi=72)
plt.show()

# 標高を緯度方向に線形補間しながら拡張する
depthsY = depths[:,0]

linearf = interpolate.interp1d(WordLats, depthsY, kind='linear', fill_value='extrapolate')
depthsExtendY = linearf(pointY)

plt.figure(figsize=(12, 9))
plt.plot(WordLats, depthsY, "o", label="original")
plt.plot(pointY, depthsExtendY, "r", label="linear")

plt.legend()
plt.savefig('2.jpg', dpi=72)
plt.show()

実験結果は以下で

経度方向に（東経110度上）の標高データをオリジナルより要素数を増やして、足りない部分は線形補間で値を補いながら結果をグラフにしてみるとこんな感じ。

青丸はオリジナルの標高・水深データで、ピクセル座標単位（元のデーターはピクセル座標ピッタリではありません）に線形補間でデータ生成してみた結果が赤線なりますが、なんかいい感じですね。しかも、オリジナルの隣接するデータからの距離も加味して補間さてれますのでなんか、使えそうですね〜〜。

なお、横軸はピクセル座標、縦軸はメートル単位の標高・水深の値で、図的には、西沙諸島の西の南シナ海からフィリピンルソン島を通り、フィリピン海、北マリアナ諸島を通り、南太平洋のウェーク島辺りまで、西から東に東経110度線上の標高・水深データーの断面になります。

Pythonで線形補間　（標高メッシュデーター経度方向に拡張してみる）

で、こちらは緯度方向に（北緯16度上）で同様の処理をしてみたもので、バイカル湖の北ぐらいから西沙諸島の西の南シナ海まで北から南に北緯16度線上の標高・水深データーの断面になります。

Pythonで線形補間　（標高メッシュデーター緯度方向に拡張してみる）

ということで、これは使えそうということで、X軸方向に線形補間で標高メッシュを拡張してから、それをさらにY軸方向に拡張してあげれば、任意のズームレベルのそれらしいメッシュデーターが作れそうということで、実際に試してみた結果がこちら。

いい感じに描けてますね。

地図の専門家からしたら、補間の仕方はまだまだでしょうが、私の個人的使用には十分で、結果的に、隣接する東西南北４点の垂直・水平方向の距離も加味しての補間になっている感じでしょうか。　今回は横・縦の順で補間しましたが、縦・横順の補間結果も生成して、２つを合算して２で割れば、より補間精度はあがるのかな？

図は、生成したズームレベル７の標高データーを使い、縦横5x5ぐらいのタイルデータを使用し地図を作成してみたもの。

こちらは、同じくズームレベル７のデータで縦横3x3タイルを使い地図にしてみたもの。

データを補間してみた時のソースは以下です

#!/usr/bin/env python
# -*- coding: utf-8 -*- 

from math import pi
from math import tanh
from math import sin
from math import asin
from numpy import arctanh

from earthquake_def import *
import MySQLdb

from netCDF4 import Dataset
import numpy as np

from scipy import interpolate



# 経緯度のピクセル座標を求める
def fromLatLngToPoint(lat,lon, z):
    L = 85.05112878
    pointX = int( (2**(z+7)) * ((lon/180) + 1) )
    pointY = int( (2**(z+7)/pi) * (-1 * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180))) )
    return pointX, pointY


# 経緯度の世界座標を求める
def fromLatLngToWordLatLon(lat, lon, z):
    L = 85.05112878
    wordLon = (2**(z+7)) * ((lon/180) + 1)
    wordLat = (2**(z+7)/pi) * (-1. * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180)))
    return wordLat, wordLon
    

# ズームレベル
z = 7

#地形データを読み込む
etopo1 = Dataset('ETOPO1_Bed_g_gmt4.grd','r')

lons   = etopo1.variables['x'][:]
lats   = etopo1.variables['y'][:]
depths = etopo1.variables['z'][:]


print lats.shape
print lons.shape
print depths.shape

# 経緯度抽出条件
condLons = (lons >= 110) * (lons <= 170)
condLats = (lats >= 16) * (lats <= 56)

# 経緯度条件を満たす配列インデックスを抽出
lonsIX = np.array(np.where(condLons)).flatten()
latsIX = np.array(np.where(condLats)).flatten()

# 経緯度条件を満たす経緯度と標高データを取り出す
lons = lons[lonsIX].flatten()
lats = lats[latsIX].flatten()
depths = depths[latsIX[0]:latsIX[-1]+1, lonsIX[0]:lonsIX[-1]+1]

# 経緯度を世界座標に変換する
WordLats = np.array([fromLatLngToWordLatLon(lat,0, z)[0] for lat in lats])
WordLons = np.array([fromLatLngToWordLatLon(0,lon, z)[1] for lon in lons])

xdiff =  int(WordLons[-1]) - int(WordLons[0])
ydiff =  int(WordLats[0]) - int(WordLats[-1])

print WordLats[0], WordLats[-1], ydiff
print WordLons[0], WordLons[-1], xdiff


# 拡張したい経度方向、緯度方向の座標を生成
# 補完後の座標がピクセル座標ピッタリになるよう間隔は1づつの整数値になるようにする
pointX = np.linspace(int(WordLons[0]), int(WordLons[-1]), xdiff+1)
pointY = np.linspace(int(WordLats[0]), int(WordLats[-1]), ydiff+1)



# 標高を経度方向に線形補間しながら拡張する
for y in range(depths.shape[0]):
    depthsX = depths[y,:]
    linearf = interpolate.interp1d(WordLons, depthsX, kind='linear', fill_value='extrapolate')
    
    # 縦に連結
    depthsLinear = [linearf(pointX)]
    if y == 0:
        depthsExtendX = depthsLinear
    else:
        depthsExtendX = np.concatenate((depthsExtendX, depthsLinear), axis = 0) #縦



# 標高を緯度方向に線形補間しながら拡張する
for x in range(depthsExtendX.shape[1]):
    depthsY = depthsExtendX[:,x]
    linearf = interpolate.interp1d(WordLats, depthsY, kind='linear', fill_value='extrapolate')
    
    # 横に連結
    depthsLinear = np.array([linearf(pointY)]).T  # 結果は1行n列からn行1列に転置しておく
    if x == 0:
        depthsExtend = depthsLinear
    else:
        depthsExtend = np.concatenate((depthsExtend, depthsLinear), axis = 1) #横


conn = MySQLdb.connect(db=e_db, host=e_db_host, port=e_db_port, user=e_db_user, passwd=e_db_passwd)
cursor = conn.cursor()

for x in range(pointX.size):
   for y in range(pointY.size):
       print x, y, depthsExtend[y, x]

       sql = '''insert IGNORE into dem_%02d set
                                     x = %d,
                                     y = %d,
                                     depth = %f
                                 ''' \
                                 % ( 
                                     z,
                                     pointX[x],
                                     pointY[y],
                                     depthsExtend[y, x]
                                  )
                                     
       print sql
       cursor.execute( sql )

cursor.close()
conn.close()

プログラム日記φ(．．)

おもにPython関係のプログラムメモ

Pythonで線形補間（欠けている標高データを補完してみる）