標高地図の作成、標高データの処理、海・陸の判断（北陸地方）

Last updated: 2022/07/29

1. 標高データ

標高データのダウンロードにはもともと専用のソフトが必要です。データのフーマットも特殊なので、その処理も専用のソフトが必要です。それぞれのソフトのインストール等は結構手間がかかるので、ここは既にダウンロードできた標高データが汎用のフーマットに変換されたあとのデータを提供します。専用のソフトが必要なく、pythonで直接処理できます（python3, Spyderの環境がおすすめです）。ただし、以下のデータは緯度35.320977〜38.199023、経度134.963775〜138.556225範囲内のデータのみです。北陸地方全体は概ね含まれています。

1.1 データのダウンロード

以下のリンクでダウンロードしてください（5.5MB）。
https://tingwu.info/pylab/data/Hokuriku_-160_-160_160_160.bz2

1.2 データの解凍

以下のコードでbz2データを解凍してください。
データ保存場所だけはダウンロードした地形データの実際の保存場所に変更する必要があります。
実行されたあと、同じフォルダーにHokuriku_-160_-160_160_160.npyが作られます。最初のbz2データを削除しても構いません。

import bz2

#データ保存場所（フールパス）
mapDataFile = 'D:/temp/Hokuriku_-160_-160_160_160.bz2'

#Decompress
tempMapFile = mapDataFile[0:-3]+'npy'
zipfile = bz2.BZ2File(mapDataFile)
data = zipfile.read()
fid = open(tempMapFile, 'wb')
fid.write(data)
fid.close()
import bz2

#データ保存場所（フールパス）
mapDataFile = 'D:/temp/Hokuriku_-160_-160_160_160.bz2'

#Decompress
tempMapFile = mapDataFile[0:-3]+'npy'
zipfile = bz2.BZ2File(mapDataFile)
data = zipfile.read()
fid = open(tempMapFile, 'wb')
fid.write(data)
fid.close()

2. 緯度経度座標の標高地図作成

以下のコードに変更する必要なところ：

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcol

#データ保存場所（フールパス）
mapDataFile = 'D:/temp/Hokuriku_-160_-160_160_160.npy' #行7

#表示する範囲
showLatiMin = 36.26 #行10 >35.320977
showLatiMax = 37.06 # <38.199023
showLongiMin = 136.26 # >134.963775
showLongiMax = 137.06 # <138.556225

#Max ranges (Do not change)
latiMin = 35.320977
latiMax = 38.199023
longiMin = 134.963775
longiMax = 138.556225

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)
X = np.linspace(longiMin, longiMax, xLen)
Y = np.linspace(latiMin, latiMax, yLen)

#Make colors to show the map
m = np.arange(0, 7.9, 0.2)
bounds = [0]
for val in m:
 	bounds.append(np.exp(val))
 	
cmap = plt.cm.gist_earth
cmaplist = [cmap(i) for i in range(cmap.N)]
del cmaplist[-5:]
del cmaplist[0:5]
cmap = matplotlib.colors.LinearSegmentedColormap.from_list('mcm',cmaplist, len(bounds)+1)

cmap_bounds = np.array(bounds)
norm = mcol.BoundaryNorm(cmap_bounds,cmap.N)

#Make figure
te = 0.03
be = 0.09
h = 1-te-be
le = 0.1
re = 0.1
hb = 0.02
w2 = 0.04
w1 = 1-le-re-hb-w2

f = plt.figure(figsize = (7.2, 6))
ax1 = f.add_axes([le, be, w1, h])
ax2 = f.add_axes([le+w1+hb, be, w2, h])
im = ax1.contourf(X, Y, data_array, bounds, cmap=cmap, norm=norm)

ax1.set_xlim([showLongiMin, showLongiMax])
ax1.set_ylim([showLatiMin, showLatiMax])

cb = plt.colorbar(im, cax=ax2)
cb.set_label('Meters above Sea Level')
ax1.set_xlabel('Longitude (Degree East)')
ax1.set_ylabel('Latitude (Degree North)')

cb.set_ticks([0, 2, 5, 10, 25, 50, 100, 250, 500, 1000, 2000])

plt.savefig('D:/temp/map.png', dpi=300)	 #行68

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcol

#データ保存場所（フールパス）
mapDataFile = 'D:/temp/Hokuriku_-160_-160_160_160.npy' #行7

#表示する範囲
showLatiMin = 36.26 #行10 >35.320977
showLatiMax = 37.06 # <38.199023
showLongiMin = 136.26 # >134.963775
showLongiMax = 137.06 # <138.556225

#Max ranges (Do not change)
latiMin = 35.320977
latiMax = 38.199023
longiMin = 134.963775
longiMax = 138.556225

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)
X = np.linspace(longiMin, longiMax, xLen)
Y = np.linspace(latiMin, latiMax, yLen)

#Make colors to show the map
m = np.arange(0, 7.9, 0.2)
bounds = [0]
for val in m:
 	bounds.append(np.exp(val))
 	
cmap = plt.cm.gist_earth
cmaplist = [cmap(i) for i in range(cmap.N)]
del cmaplist[-5:]
del cmaplist[0:5]
cmap = matplotlib.colors.LinearSegmentedColormap.from_list('mcm',cmaplist, len(bounds)+1)

cmap_bounds = np.array(bounds)
norm = mcol.BoundaryNorm(cmap_bounds,cmap.N)

#Make figure
te = 0.03
be = 0.09
h = 1-te-be
le = 0.1
re = 0.1
hb = 0.02
w2 = 0.04
w1 = 1-le-re-hb-w2

f = plt.figure(figsize = (7.2, 6))
ax1 = f.add_axes([le, be, w1, h])
ax2 = f.add_axes([le+w1+hb, be, w2, h])
im = ax1.contourf(X, Y, data_array, bounds, cmap=cmap, norm=norm)

ax1.set_xlim([showLongiMin, showLongiMax])
ax1.set_ylim([showLatiMin, showLatiMax])

cb = plt.colorbar(im, cax=ax2)
cb.set_label('Meters above Sea Level')
ax1.set_xlabel('Longitude (Degree East)')
ax1.set_ylabel('Latitude (Degree North)')

cb.set_ticks([0, 2, 5, 10, 25, 50, 100, 250, 500, 1000, 2000])

plt.savefig('D:/temp/map.png', dpi=300)	 #行68

3. FALMA直角座標の標高地図作成

2のコードと同じように、以下のところが変更する必要があります：

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcol

#データ保存場所（フールパス）
mapDataFile = 'D:/temp/Hokuriku_-160_-160_160_160.npy' #行7

#表示する範囲 (km)
showXMin = -50 #行10 >-160
showXMax = 50 # <160
showYMin = -50 # >-160
showYMax = 50 # <160

#Max ranges (Do not change)
xMin = -160
xMax = 160
yMin = -160
yMax = 160

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)
X = np.linspace(xMin, xMax, xLen)
Y = np.linspace(yMin, yMax, yLen)

#Make colors to show the map
m = np.arange(0, 7.9, 0.2)
bounds = [0]
for val in m:
 	bounds.append(np.exp(val))
 	
cmap = plt.cm.gist_earth
cmaplist = [cmap(i) for i in range(cmap.N)]
del cmaplist[-5:]
del cmaplist[0:5]
cmap = matplotlib.colors.LinearSegmentedColormap.from_list('mcm',cmaplist, len(bounds)+1)

cmap_bounds = np.array(bounds)
norm = mcol.BoundaryNorm(cmap_bounds,cmap.N)

#Make figure
te = 0.03
be = 0.09
h = 1-te-be
le = 0.1
re = 0.1
hb = 0.02
w2 = 0.04
w1 = 1-le-re-hb-w2

f = plt.figure(figsize = (7.2, 6))
ax1 = f.add_axes([le, be, w1, h])
ax2 = f.add_axes([le+w1+hb, be, w2, h])
im = ax1.contourf(X, Y, data_array, bounds, cmap=cmap, norm=norm)

ax1.set_xlim([showXMin, showXMax])
ax1.set_ylim([showYMin, showYMax])

cb = plt.colorbar(im, cax=ax2)
cb.set_label('Meters above Sea Level')
ax1.set_xlabel('West-East (km)')
ax1.set_ylabel('South-North (km)')

cb.set_ticks([0, 2, 5, 10, 25, 50, 100, 250, 500, 1000, 2000])

plt.savefig('D:/temp/map.png', dpi=300)	#行68

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as mcol

#データ保存場所（フールパス）
mapDataFile = 'D:/temp/Hokuriku_-160_-160_160_160.npy' #行7

#表示する範囲 (km)
showXMin = -50 #行10 >-160
showXMax = 50 # <160
showYMin = -50 # >-160
showYMax = 50 # <160

#Max ranges (Do not change)
xMin = -160
xMax = 160
yMin = -160
yMax = 160

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)
X = np.linspace(xMin, xMax, xLen)
Y = np.linspace(yMin, yMax, yLen)

#Make colors to show the map
m = np.arange(0, 7.9, 0.2)
bounds = [0]
for val in m:
 	bounds.append(np.exp(val))
 	
cmap = plt.cm.gist_earth
cmaplist = [cmap(i) for i in range(cmap.N)]
del cmaplist[-5:]
del cmaplist[0:5]
cmap = matplotlib.colors.LinearSegmentedColormap.from_list('mcm',cmaplist, len(bounds)+1)

cmap_bounds = np.array(bounds)
norm = mcol.BoundaryNorm(cmap_bounds,cmap.N)

#Make figure
te = 0.03
be = 0.09
h = 1-te-be
le = 0.1
re = 0.1
hb = 0.02
w2 = 0.04
w1 = 1-le-re-hb-w2

f = plt.figure(figsize = (7.2, 6))
ax1 = f.add_axes([le, be, w1, h])
ax2 = f.add_axes([le+w1+hb, be, w2, h])
im = ax1.contourf(X, Y, data_array, bounds, cmap=cmap, norm=norm)

ax1.set_xlim([showXMin, showXMax])
ax1.set_ylim([showYMin, showYMax])

cb = plt.colorbar(im, cax=ax2)
cb.set_label('Meters above Sea Level')
ax1.set_xlabel('West-East (km)')
ax1.set_ylabel('South-North (km)')

cb.set_ticks([0, 2, 5, 10, 25, 50, 100, 250, 500, 1000, 2000])

plt.savefig('D:/temp/map.png', dpi=300)	#行68

地図の上に、観測サイトの場所や雷の標定点をプロットすることができます。

3. 緯度経度、FALMA座標から標高の算出

標高のデータがあるので、緯度経度やFALMA座標の情報があれば、標高をすぐ算出できます。

3.1 緯度経度から標高を算出

import numpy as np

#データ保存場所（フールパス）
mapDataFile = '/home/wu/temp/Hokuriku_-160_-160_160_160.npy'

#Max ranges (Do not change)
latiMin = 35.320977
latiMax = 38.199023
longiMin = 134.963775
longiMax = 138.556225

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)

def getElevation(latiVal, longiVal):
    latiIndex = int(round((latiVal-latiMin)/(latiMax-latiMin)*yLen))
    longiIndex = int(round((longiVal-longiMin)/(longiMax-longiMin)*xLen))
    return data_array[latiIndex, longiIndex]

elevationVal = getElevation(36.1424992, 136.7684124) #行22
print(elevationVal)
import numpy as np

#データ保存場所（フールパス）
mapDataFile = '/home/wu/temp/Hokuriku_-160_-160_160_160.npy'

#Max ranges (Do not change)
latiMin = 35.320977
latiMax = 38.199023
longiMin = 134.963775
longiMax = 138.556225

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)

def getElevation(latiVal, longiVal):
    latiIndex = int(round((latiVal-latiMin)/(latiMax-latiMin)*yLen))
    longiIndex = int(round((longiVal-longiMin)/(longiMax-longiMin)*xLen))
    return data_array[latiIndex, longiIndex]

elevationVal = getElevation(36.1424992, 136.7684124) #行22
print(elevationVal)

行22に緯度経度の値を入れて、実行すれば、その緯度経度のところの標高がプリントされる。

3.2 FALMAの座標から標高を算出

import numpy as np

#データ保存場所（フールパス）
mapDataFile = '/home/wu/temp/Hokuriku_-160_-160_160_160.npy' 

#Max ranges (Do not change)
xMin = -160
xMax = 160
yMin = -160
yMax = 160

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)

def getElevation(xVal, yVal):
    xIndex = int(round((xVal-xMin)/(xMax-xMin)*xLen))
    yIndex = int(round((yVal-yMin)/(yMax-yMin)*yLen))
    return data_array[yIndex, xIndex]

elevationVal = getElevation(10, -20) #行22
print(elevationVal)
import numpy as np

#データ保存場所（フールパス）
mapDataFile = '/home/wu/temp/Hokuriku_-160_-160_160_160.npy' 

#Max ranges (Do not change)
xMin = -160
xMax = 160
yMin = -160
yMax = 160

#Load data
data_array = np.load(mapDataFile)

[yLen, xLen] = np.shape(data_array)

def getElevation(xVal, yVal):
    xIndex = int(round((xVal-xMin)/(xMax-xMin)*xLen))
    yIndex = int(round((yVal-yMin)/(yMax-yMin)*yLen))
    return data_array[yIndex, xIndex]

elevationVal = getElevation(10, -20) #行22
print(elevationVal)

行22にFALMA座標(km)の値を入れて、実行すれば、その座標の標高がプリントされる。

4. 海・陸の判断

3のコードを利用して、もしelevationValの結果はnanであれば、その緯度経度・座標のところは海と判断できます。
例えば、3.1のコードを利用して、(36.7777058, 136.7114041)の標高を計算すれば、結果はnanなので、この緯度経度のところは海ということが判断できます。
ただし、海岸線付近数百メートル以内は正確に判断できない場合もあります。いろんなところの標高を計算して、GoogleMapと比較してみてください。

Back to Python関連資料