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

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

Pythonでコピーフリーの海底地形き日本地図を作ってみる

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

ETOPO1とOpenStreetMapで作るコピーフリーの海底地形図付き日本地図

ETOPO1とOpenStreetMapで作るコピーフリーの海底地形図付き日本地図

 

作り方

 

ETOPO1グローバル地形データセットの説明は以前に書いた以下の記事や

 

memomemokun.hateblo.jp

 

タイル形式のOpenStreetMapの地図データの地図タイルの考え方については以下の記事などを参照いただくとして

 

memomemokun.hateblo.jp

 

OpenStreetMapの地図タイル画像にETOPO1グローバル地形データーを乗せ易いように、ETOPO1の地形(標高・水深)データをピクセル座標のメッシュ単位の標高・水深データーに変換してMySQLに一旦格納。

 

#!/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

import MySQLdb
from io import StringIO

from netCDF4 import Dataset
import pandas as pd
import numpy as np

e_db_host = "localhost"
e_db_port = 3306
e_db = "earthquake"
e_db_user = "user"
e_db_passwd = "passwd"


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

# ズームレベルに合わせ以下は変更
z = 6

# ETOPO1の地形データを読み込み、ピクセル座標を付加し変数 maptxt に格納
etopo1 = Dataset('ETOPO1_Bed_g_gmt4.grd','r')

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


lats = np.array(lats)
lons = np.array(lons)
depths = np.array(depths)

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

f = StringIO()
for x in range(depths.shape[1]):
    for y in range(depths.shape[0]):

        # データーは日本周辺に絞る
        if lons[x] < 110 or lons[x] > 170: continue
        if lats[y] < 16 or lats[y] > 56: continue
        pointX, pointY = fromLatLngToPoint(lats[y], lons[x], z)
    
        #print x, y, lons[x], lats[y], depths[y,x], pointX, pointY 
        txt = u'{0:.0f},{1:.0f},{2:.12f},{3:.12f},{4:.0f}'.format(pointX, pointY, lons[x], lats[y], depths[y,x])
        f.write(txt + "\n")
    
maptxt = f.getvalue()
f.close()


#変数 maptxtから地形データを読み込む
jodc = pd.read_csv(StringIO(maptxt), header=None, 
                         names=('x', 'y', 'lat', 'lon', 'depth'), 
                         dtype={'col1':'i1', 'col2':'i1', 'col3':'f8', 'col3':'f8', 'col5':'i1'})


#深度データ最大、最小値などを求めておく
xmin = jodc['x'].min()
xmax = jodc['x'].max()
xdiff = jodc['x'].max()-jodc['x'].min()
ymin = jodc['y'].min()
ymax = jodc['y'].max()
ydiff = jodc['y'].max()-jodc['y'].min()
depth = jodc['depth'].min()

print xmin, xmax, xdiff
print ymin, ymax, ydiff

#深度データの点群データを格子データに変換する
Density = np.zeros( (ydiff+1, xdiff+1) ) #格子内の点の数
Intensity = np.zeros( (ydiff+1, xdiff+1) ) #Z値の総和
for index, row in jodc.iterrows():
    #print row.x-xmin, row.y-ymin, row.depth
    x = int(row.x-xmin)
    y = int(row.y-ymin)
    Density[y,x] = Density[y,x] + 1
    Intensity[y,x] = Intensity[y,x] + row.depth

#格子ごとのZ値の総和 / 格子内の点の数 として各格子単位の深度データを求める
Z=Intensity/Density

#欠損値は取り敢えず0にする
Z[np.isnan(Z)] = 0.0

print Z.shape

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(0,xdiff+1):
   for y in range(0,ydiff+1):
       print x, y, Z[y, x]

       sql = '''insert IGNORE into dem_%02d set
                                     x = %d,
                                     y = %d,
                                     depth = %f
                                 ''' \
                                 % ( 
                                     z,
                                     xmin+x,
                                     ymin+y,
                                     Z[y, x]
                                  )
                                     
       print sql
       cursor.execute( sql )

cursor.close()
conn.close()
    

 

 

で、

 

後は、OpenStreetMap国土地理院の標高データーを合体して地形の陰影付き地図を作成した時の要領で、MySQLに格納したピクセル座標単位のETOPO1の地形(標高・水深)データをOpenStreetMapの地図画像を合成してあげると

 

memomemokun.hateblo.jp

 

 

こんな、感じの海底地形図付き日本地図が出来上がります。

画像は、北緯36度、東経141度の地点を中心に、北緯26度、東経130度から北緯46度、東経152度の範囲をズームレベル5のメッシュ精度で縦横5X5の地図タイル合成して日本付近を描いたもの。

 

南はフィリピン、台湾付近から、中国、朝鮮半島、シベリア、カムチャッカから伊豆・小笠原からマリアナ海溝まで、、オホーツク海、太平洋、フリピン海の海溝や海山などの海底地形も描けてますね。

 

地図画像にはOpenStreetMapの著作権とライセンスの規定に従い、生成した地図画像には「© OpenStreetMap contributors」のクレジットも付加。

 

(以下は1,728x1,728サイズの拡大地図ので、興味のある方は拡大してみてください。)

ETOPO1とOpenStreetMapで作るコピーフリーの海底地形図付き日本地図

ETOPO1とOpenStreetMapで作るコピーフリーの海底地形図付き日本地図

 

 

なお、ETOPO1の地形データは全地球で2億個程もあるメッシュデーターですが、地球は広く、南極・北極の極域も含め全地球の海陸あわせた1分=約1.8kmの解像度の標高・水深のグリッドデータですので、拡大すると、ズームレベル6ぐらいが精一杯でしょうか。

 

ETOPO1の地形データを線形補間か2 次スプライン補間でもかけてデータが無いグリットはそれらしいデータを作ってやれば2ズームぐらいは引き伸ばせるかもしれませんね。

 

というか、実際に線形補間で足りない部分はそれらしい値でメッシュを補ってズーム7で陰影起伏を描いてみたのが以下で、やってみると、ズーム12ぐらいに引き伸ばしても、それらしい、陰影起伏は描けるみたいですね。

 

memomemokun.hateblo.jp

 

なお、今回標高データを使用したのは、陰影起伏図と海底の水深の色分けのために使用しただけですので、最初からズームレベル6の陰影起伏図と色分けした水深画像を作成しておいて、画像を3倍ぐらいに拡大しても、それ以上の解像度は望めませんが地形の雰囲気とかは伝わると思うので、最初から、ズーム9までの陰影起伏図と色分けした水深画像の画像タイルを拡大画像も利用し作成しておけば、あとは、単に画像の合成に帰着しちゃいますね。

 

北緯36度、東経141度の地点を中心に、縦横2x2タイルで作成してみたもの。

f:id:memomemokun:20190112151211j:plain

 

 

さらに細かい精度で描きたいときは

 

地上部は、スペースシャトルの標高データSRTMや地球観測衛星「だいち」の標高モデル

http://memomemokun.hateblo.jp/entry/2018/12/09/231706

http://memomemokun.hateblo.jp/entry/2019/01/15/204737

 

海洋部は、海上保安庁日本海洋データセンター(JODC)の日本周辺の500mグリッド

http://www.jodc.go.jp/data_set/jodc/jegg_intro_j.html

http://memomemokun.hateblo.jp/entry/2018/11/12/190414

 

などを合わせて行くと、より細かく描けるでしょうかね。著作権の確認は必要ですが。

 

ということで、上記を描画した際のソースは以下

 

以下ソースの説明に、以下記事も参照いただくとやってることがわかるかも。

 

memomemokun.hateblo.jp

 

memomemokun.hateblo.jp

 

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
   
import os
import requests
from StringIO import StringIO
import string
import pandas as pd
import numpy as np
import hjson

import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib import cm

from scipy import misc

from mpl_toolkits.basemap import Basemap
from osgeo import gdal

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

from pyproj import Geod

import datetime
from dateutil.relativedelta import relativedelta
import matplotlib.font_manager as fm


import MySQLdb

e_db_host = "localhost"
e_db_port = 3306
e_db = "earthquake"
e_db_user = "user"
e_db_passwd = "passwd"

prop_j = fm.FontProperties(fname='/system/Library/Fonts/ヒラギノ丸ゴ ProN W4.ttc')

# 1タイルのピクセル数
pix = 256

# 描画に使用する縦横の最小タイル数
Tiles = 5

g = Geod(ellps='WGS84')


# ピクセル座標を経緯度に変換する
def fromPointToLatLng(pixelLat,pixelLon, z):
    L = 85.05112878
    lon = 180 * ((pixelLon / 2.0**(z+7) ) - 1)
    lat = 180/pi * (asin(tanh(-pi/2**(z+7)*pixelLat + arctanh(sin(pi/180*L)))))
    return lat, lon


# 経緯度をピクセル座標に変換する
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 fromLatLngToTile(lat,lon, z):
    pointX, pointY = fromLatLngToPoint(lat,lon, z)
    return pointX/pix, pointY/pix    


# 指定経緯度範囲が異なるタイルにまたがる最小のズームレベルとそのタイル座標を返す
def getTile(Tiles, minLat, minLon, maxLat, maxLon):
    
    for z in range(0,7):
        tileX1, tileY1 = fromLatLngToTile(minLat, minLon, z)
        tileX2, tileY2 = fromLatLngToTile(maxLat, maxLon, z)
        
        if tileX2 - tileX1 >= Tiles - 1  and tileY1 - tileY2 >= Tiles - 1: break

    return z, tileX1, tileY1, tileX2, tileY2 



# 指定緯度帯のZ係数を返す
def Z_factor(lat):

    zfactors = {
             0 : 0.00000898,
            10 : 0.00000912,
            20 : 0.00000956,
            30 : 0.00001036,
            40 : 0.00001171,
            50 : 0.00001395,
            60 : 0.00001792,
            70 : 0.00002619,
            80 : 0.00005156
        }
    
    return zfactors[np.trunc(lat/10.0)*10.0]



# 地図タイルを読み込む
def readTile(tiles, tilename, x, y, z):

    # タイル用ディレクトリを作成する
    if not os.path.isdir('tiles'): os.mkdir('tiles')
    
    tiledir = 'tiles/%s' % tilename
    if not os.path.isdir(tiledir): os.mkdir(tiledir)
    
    tiledir = 'tiles/%s/%d' % (tilename, z)
    if not os.path.isdir(tiledir): os.mkdir(tiledir)
    
    tiledir = 'tiles/%s/%d/%d' % (tilename, z, x)
    if not os.path.isdir(tiledir): os.mkdir(tiledir)
    
    #tiledir = 'tiles/%s/%d/%d/%d' % (tilename, z, x, y)
    #if not os.path.isdir(tiledir): os.mkdir(tiledir)
    
    # タイルファイル名
    extension =  os.path.splitext(tiles['tiles'][tilename])[1][1:]
    tilef = 'tiles/%s/%d/%d/%d.%s' % (tilename, z, x, y, extension)

    # 地図タイルがローカルに存在しないときはネットからダウンロード
    if not os.path.isfile(tilef):

        url = tiles['tiles'][tilename].format(z=z,x=x,y=y)
        print url
        response = requests.get(url)
        if response.status_code == 404:

            if extension == 'txt':
            
                #標高データが無い区画は'e'で書き出す
                ndarr = np.array(['e']*256*256).reshape( (256,256) )
                np.savetxt(tilef, ndarr, delimiter=',', fmt='%s')

            else:

                #地図画像データが無い区画は白塗画像を書き出す
                img=np.ones((256, 256, 3), dtype=object)
                img *= 255
                img = np.uint8(img)
                misc.imsave(tilef, img)

        else:

            if extension == 'txt':

                # テキストダウンロード
                f=open(tilef, mode='w')
                f.write(response.content)
                f.close()

            else:

                # 画像ダウンロード
                img = misc.imread(StringIO(response.content))
                misc.imsave(tilef, img)


    # 地図画像データを読み込む
    f=open(tilef, mode='rb')
    tile = f.read()
    f.close()
    
    return tile



# 地図画像データを読み込む
def load_imgColors(tiles, tilename, z, x1, y1, x2, y2):

    for x in range(x1, x2+1):
        for y in range(y2, y1+1):
            
            #地図画像データを読み込む
            img = readTile(tiles, tilename, x, y, z)
            img = misc.imread(StringIO(img))
                 
            # 画像タイルを縦に連結
            if y == y2:
                 im_v = img
            else:
                 #im_v = cv2.vconcat([im_v, img])
                 #im_v = np.append(im_v,img,0)
                 im_v = np.concatenate((im_v, img), axis = 0) #縦
                
        # 画像タイルを横に連結
        if x == x1:
            im = im_v
        else:
            #im = cv2.hconcat([im, im_v])
            #im = np.append(im,im_v,1)
            im = np.concatenate((im,im_v), axis = 1) #横
            
    return im



# 地形データを読み込む
def load_gis(z, x1, y1, x2, y2):
    
    pixX1 = x1 * pix
    pixY1 = y2 * pix
    pixX2 = x1 * pix + 255 * (x2-x1+1)
    pixY2 = y2 * pix + 255 * (y1-y2+1)
    
    print pixX1, pixX2
    print pixX2, pixY2

    # 水深メッシュデータを読み込む
    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()
    sql = "select * from dem_%02d where x>=%d and x<=%d and y>=%d and y<=%d" % (z, pixX1, pixX2, pixY1, pixY2)
    #print sql
    cnt = cursor.execute( sql )
    rec = cursor.fetchall()
    cursor.close()
    conn.close()
    
    gis = np.zeros( (256 * (y1-y2+1), 256 * (x2-x1+1)) )
    for item in rec:
        gis[item[1]-pixY1, item[0]-pixX1] = item[2]
        
    
    return gis



# 標高タイルから陰影起伏 (Hillshade) を計算する
def hillshade(Z, center_Lat, azimuth_deg, altitude_deg):

    # 通常の陰影起伏図は、光源方位 315 度、光源方位 45 度で描画
    #
    # 陰影起伏 = 255.0 * ((cos(Zenith_rad) * cos(Slope_rad)) +
    #                 (sin(Zenith_rad) * sin(Slope_rad) * cos(Azimuth_rad - Aspect_rad)))    


    # 勾配(傾き)を求める
    Zy, Zx = np.gradient(Z)

    # 傾斜角
    Slope_rad = np.arctan(30 * Z_factor(center_Lat) * np.sqrt(Zx*Zx + Zy*Zy))

    # 傾斜方向
    Aspect_rad = np.arctan2(-Zx, Zy)

    # 太陽の方位角をラジアンに変換
    Azimuth_rad = np.radians(360.0 - azimuth_deg + 90)

    # 太陽の天頂角をラジアンに変換
    Zenith_rad = np.radians(90 - altitude_deg)

    # 陰影起伏
    shaded = np.cos(Zenith_rad) * np.cos(Slope_rad) + np.sin(Zenith_rad) * np.sin(Slope_rad) * np.cos(Azimuth_rad - Aspect_rad)
    return 255 * shaded



# 画像を平坦化する
def imgFlattening(img):
    
    # 画像のヒストグラムを平坦化する
    # ヒストグラムを取得
    hist, bins = np.histogram( img.flatten(), bins=256 )  

    # 累積分布を取る  
    cdf = hist.cumsum()

    # 0〜255の範囲に正規化
    cdf = 255 * cdf / cdf[-1]

    # ヒストグラムが平坦化するよう線形補間
    img2 = np.interp( img.flatten(), bins[:-1], cdf)

    return img2.reshape( img.shape )


def makeMap(center_Lat, center_Lon, minLat, minLon, maxLat, maxLon):

    # 指定経緯度範囲の最適ズームとタイル座標範囲を求める
    z, x1, y1, x2, y2  = getTile(Tiles, minLat, minLon, maxLat, maxLon)
    print z, x1, y1, x2, y2

    
    # タイルURL
    f=open('tiles.hjson')
    tiles = hjson.load(f)
    f.close()


    # OpenStreetMapのタイル画像を読み込み連結して1枚の画像として返す
    imgColors = load_imgColors(tiles, 'OSMA', z, x1, y1, x2, y2) 

    # 標高タイルを読み込み連結して1枚の標高タイルとして返す
    Z = load_gis(z, x1, y1, x2, y2)


    # 標高タイルから陰影起伏を計算する
    Zshade = hillshade(Z, center_Lat, 315, 60)
        
    # 陰影起伏を0..1に正規化
    Zshade_norm = (Zshade-Zshade.min())/(Zshade.max()-Zshade.min())

    # 陰影起伏をグレイスケール画像にする
    projectionIMG = cm.binary(Zshade_norm)*255


    # 透過情報はカット
    projectionIMG = projectionIMG[:,:,:3]
    #projectionIMG = np.where(projectionIMG == 125, 255, np.uint8(projectionIMG))

    # 陰影起伏図の画像を平坦化して明度のバランスを整える
    projectionIMG = imgFlattening(projectionIMG)

    
    # なんとなくガンマ補正 
    gamma = 1.6
    imax = projectionIMG.max() 
    projectionIMG = imax * (projectionIMG / imax)**(1/gamma)


    # 背景地図を一旦明度調整
    imgColors = imgFlattening(imgColors)

    # 地図画像と射影印影図を足して2で割りゃ画像が合成される
    imgColors = (imgColors + projectionIMG) / 2.0
 

    imgColors = np.uint8(imgColors)

    # 指定経緯度のピクセル座標を求める
    pointX1, pointY1 = fromLatLngToPoint(minLat, minLon, z)
    pointX2, pointY2 = fromLatLngToPoint(maxLat, maxLon, z)


    # 経緯度範囲に地図タイルをカット
    imgColors = imgColors[pointY2-y2*pix:pointY2-y2*pix+pointY1-pointY2, pointX1-x1*pix:pointX1-x1*pix+pointX2-pointX1, :]
    projectionIMG = projectionIMG[pointY2-y2*pix:pointY2-y2*pix+pointY1-pointY2, pointX1-x1*pix:pointX1-x1*pix+pointX2-pointX1, :]
    Z = Z[pointY2-y2*pix:pointY2-y2*pix+pointY1-pointY2, pointX1-x1*pix:pointX1-x1*pix+pointX2-pointX1]


    # 地図を作成する
    # 3857 球面メルカトル
    # 3395 回転楕円体メルカトル
    fig = plt.figure(figsize=(12,12))
    plt.subplots_adjust(left=0.0, right=1.0, bottom=0.04, top=0.96)
    map = Basemap(epsg=3395, 
            lat_ts=center_Lat, resolution='h',
            llcrnrlon=minLon, llcrnrlat=minLat, 
            urcrnrlon=maxLon, urcrnrlat=maxLat)
            
           

    # 背景地図画像を描画
    mx0, my0 = map(minLon, minLat)
    mx1, my1 = map(maxLon, maxLat)

    extent = (mx0, mx1, my0, my1)
    plt.imshow(imgColors, vmin = 0, vmax = 255, extent = extent, alpha=1.0)


    # 水深・標高用カラーマップ作成
    cmap = colors.ListedColormap(['blue', 'royalblue', 'steelblue', 'cornflowerblue', 'lightblue', 'paleturquoise', 'lightcyan', 'azure'])
    cmap.set_over('honeydew')
    #cmap.set_over('#c0d0c0')
    cmap.set_under('darkblue')
    bounds = [-4000, -3000, -2000, -1000, -500, -200, -100, -50, 0]
    norm = colors.BoundaryNorm(bounds, cmap.N)

            
    #X,Y軸のグリッドを生成
    x = np.linspace(0, map.urcrnrx, Z.shape[1])
    y = np.linspace(map.urcrnry, 0, Z.shape[0])
    X, Y = np.meshgrid(x, y)

    # 水深・標高に応じて背景に色を付ける
    surf = map.pcolormesh(X, Y, Z, norm=norm, alpha=0.07, cmap=cmap)

    # コピーライト表示
    x, y = map(minLon, minLat)
    plt.text(x+10000, y+30000, u'Maps © 気になるメモメモφ(..), 地図データ © OpenStreetMap contributors', fontsize=18, color='k', alpha=1.00, fontproperties=prop_j)


    # 経緯度線を描画
    map.drawparallels(np.arange(int(minLat), int(maxLat)+1, coordinateLineStep[z]), labels = [1,0,0,0], fontsize=18)
    map.drawmeridians(np.arange(int(minLon), int(maxLon)+1, coordinateLineStep[z]), labels = [0,0,0,1], fontsize=18)


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


# ズームレベルに応じた経緯度線のステップ
coordinateLineStep = {
    3 : 20.0,
    4 : 15.0,
    5 : 10.0,
    6 : 5.0
}


if __name__ == '__main__':


    center_Lat = 36.0
    center_Lon = 141.0


    gap = 20.0
    minLat = center_Lat - gap
    minLon = center_Lon - gap * 1.1
    maxLat = center_Lat + gap
    maxLon = center_Lon + gap * 1.1

    print "北緯%f 東経%f - 北緯%f 東経%f" % (minLat, minLon, maxLat, maxLon)

    makeMap(center_Lat, center_Lon, minLat, minLon, maxLat, maxLon)