Python AI 画像¶
SIGNATE社とウエザーニュース社の気象衛星画像を用いた雲画像予測コンペティションのチュートリアルを実施してみました。¶
(HP上のコードをなぞってみただけですが） https://signate.jp/competitions/169/tutorials/12
地図描画は、basemapではなく、cartopyを使用しました。
Anacondaの場合、conda install numpy、conda install matplotlib、conda install pandas、conda install keras-gpu等でインストールしてください。
condaでインストールできない場合は、pip install ternsorflow kerasのようにpipでインストールしてください。
In [1]:
import numpy as np
import cv2
import os
import gzip
import shutil # ファイル操作
from datetime import datetime as dt
from datetime import timedelta
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.ticker as tick
import warnings
warnings.filterwarnings('ignore', category=matplotlib.MatplotlibDeprecationWarning)
１．データの読み込み
In [2]:
dir = "./signate/sat_image_2017_01/"
file1 = dir + "train/sat/2017-01-01/2017-01-01-00-00.fv.png"
file4 = dir + "train/sat/2017-04-01/2017-04-01-00-00.fv.png"
file6 = dir + "train/sat/2017-06-30/2017-06-30-00-00.fv.png"
img1 = cv2.imread(file1, 0) # 2017年1月の衛星画像
img4 = cv2.imread(file4, 0) # 2017年4月の衛星画像
img6 = cv2.imread(file6, 0) # 2017年6月の衛星画像
print( "image shape : ", img1.shape )
print( "data type : ", img1.dtype )
fig, axes = plt.subplots(1, 3, figsize = (9, 6))
axes[0].imshow(img1, cmap = 'gray')
axes[0].set_title("SatelliteImage on 2017/01/01 00:00UT", fontsize = 9)
axes[1].imshow(img4, cmap = 'gray')
axes[1].set_title("SatelliteImage on 2017/04/01 00:00UT", fontsize = 9)
axes[2].imshow(img6, cmap = 'gray')
axes[2].set_title("SatelliteImage on 2017/06/30 00:00UT", fontsize = 9)
plt.show()
ヒストグラム表示
In [3]:
img01 = []
img04 = []
img06 = []
date01 = dt(2017, 1 , 1, 0, 0, 0) # 2017年1月分 年、月、日、時、分、秒
date04 = dt(2017, 4 , 1, 0, 0, 0) # 2017年4月分
date06 = dt(2017, 6 , 1, 0, 0, 0) # 2017年6月分
# 1ヶ月分の衛星画像を読み込む
for i in range(31*24):
file01 = dir + "train/sat/{year}-{month:02}-{day:02}/{year}-{month:02}-{day:02}-{hour:02}-00.fv.png".\
format(year = date01.year, month = date01.month, day = date01.day, hour = date01.hour)
if os.path.exists(file01): # ファイルの存在確認
img = cv2.imread(file01, 0)
img01.append(img) # 1月分
file04 = dir + "train/sat/{year}-{month:02}-{day:02}/{year}-{month:02}-{day:02}-{hour:02}-00.fv.png".\
format(year = date04.year, month = date04.month, day = date04.day, hour = date04.hour)
if os.path.exists(file04):
img = cv2.imread(file04, 0)
img04.append(img) # 4月分
file06 = dir + "train/sat/{year}-{month:02}-{day:02}/{year}-{month:02}-{day:02}-{hour:02}-00.fv.png".\
format(year = date06.year, month = date06.month, day = date06.day, hour = date06.hour)
if os.path.exists(file06):
img = cv2.imread(file06, 0)
img06.append(img) # 6月分
date01 = date01 + timedelta(hours = 1) # 日付の加算・
date04 = date04 + timedelta(hours = 1)
date06 = date06 + timedelta(hours = 1)
print( "image shape : ", img1.shape )
# ヒストグラム作成
all_img01 = np.concatenate(img01).flatten()
all_img04 = np.concatenate(img04).flatten()
all_img06 = np.concatenate(img06).flatten()
fig, ax = plt.subplots(1, 1, figsize = (8, 8))
ax.hist(all_img01, bins = np.arange(256 + 1), alpha = 0.3, color = "r", label = "January, 2017")
ax.hist(all_img04, bins = np.arange(256 + 1), alpha = 0.3, color = "g", label = "April, 2017")
ax.hist(all_img06, bins = np.arange(256 + 1), alpha = 0.3, color = "b", label = "June, 2017")
ax.set_title("Histogram of Satellite Images")
ax.set_xlabel("Pixel")
ax.set_ylim(0, 6000000)
plt.legend()
plt.show()
日変化画像
In [4]:
filemorning = dir + "train/sat/2017-01-31/2017-01-31-21-00.fv.png"
filenoon = dir + "train/sat/2017-01-31/2017-01-31-03-00.fv.png"
fileevening = dir + "train/sat/2017-01-31/2017-01-31-09-00.fv.png"
imgmorning = cv2.imread(filemorning, 0) # 朝の画像
imgnoon = cv2.imread(filenoon, 0) # 昼の画像
imgevening = cv2.imread(fileevening, 0) # 夕の画像
fig, axes = plt.subplots(1, 3, figsize = (9, 6))
axes[0].imshow(imgmorning, cmap = 'gray')
axes[0].set_title("2017/01/31 21:00UT", fontsize = 9)
axes[1].imshow(imgnoon, cmap = 'gray')
axes[1].set_title("2017/01/31 03:00UT", fontsize = 9)
axes[2].imshow(imgevening, cmap = 'gray')
axes[2].set_title("2017/01/31 06:00UT", fontsize = 9)
plt.show()
In [5]:
# 気象マップ表示関数
def draw_weather_map(data = None, u = None, v = None, draw = "shaded", levels = None):
# u wind 東西風
# v wind 南北風
# 緯度経度で範囲を指定する
north = 48.
south = 20.
east = 150.
west = 118.
grid_size = [168, 128] # y, x pixels
fig = plt.figure(figsize = (9, 9))
ax = fig.add_subplot(1, 1 ,1, projection = ccrs.PlateCarree()) # 正距円筒図法、標準緯線0度 ccrs.Orthographic正距円筒
# 地図の表示
ax.set_extent([west, east, south, north], ccrs.PlateCarree()) # 正距円筒図法
# 海岸線を引く
ax.coastlines(resolution = '10m', lw = 1.0, facecolor = 'lightgreen') # 解像度 50m 0.5 線が細い
# # 陸地をカラー表示
# ax.add_feature(cfeature.LAND, facecolor = 'lightgreen') # 110m resolution
land_10m = cfeature.NaturalEarthFeature('physical',
'land',
'10m',
edgecolor = 'gray',
facecolor = cfeature.COLORS['land'])
ax.add_feature(land_10m) # 10m resolution
# # グリッドを引く
# ax.gridlines(linestyle = '-', color = 'grey')
#グリッドと軸目盛を描く緯度経度を設定するための配列
dlon, dlat = 8, 7
xticks = np.arange(west, east + 0.1, dlon)
yticks = np.arange(south, north + 0.1, dlat) #目盛を描く緯度経度の値を設定
ax.set_xticks(xticks, crs = ccrs.PlateCarree())
ax.set_yticks(yticks, crs = ccrs.PlateCarree())
#目盛の表示形式を度数表記にする
latfmt = LatitudeFormatter()
lonfmt = LongitudeFormatter(zero_direction_label = True)
ax.xaxis.set_major_formatter(lonfmt)
ax.yaxis.set_major_formatter(latfmt)
#目盛のサイズを指定
ax.axes.tick_params(labelsize = 12)
# グリッドを引く
# ax.gridlines(linestyle = '-', color = 'grey')
lons = np.linspace(west, east, grid_size[1])
lats = np.linspace(south, north, grid_size[0])
gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=False,linewidth=1,alpha=0.8,linestyle='-',color='grey')
gl.xlocator = tick.FixedLocator(xticks)
gl.ylocator = tick.FixedLocator(yticks)
if draw == "shaded": # 領域
cs = ax.contourf(lons, lats, data[::-1, :], levels = levels)
fig.colorbar(cs, ax = ax, shrink = 0.7) # カラーバーのサイズを小さくした
elif draw == "contour": # 等高線
cs = ax.contour(lons, lats, data[::-1, :], levels = levels, transform = ccrs.PlateCarree())
ax.clabel(cs, fmt = '%d', fontsize = 14)
elif draw == "barbs": # 羽枝
cs = ax.contour(lons, lats, data[::-1, :], levels = levels, linewidths = 3)
ax.clabel(cs, fmt = '%d', fontsize = 14)
ax.barbs(lons, lats, # x-positions, y-positions
u, v, # u-direction, v-direction
length = 7,
sizes = dict(emptybarb = 0.25, spacing = 0.2, height = 0.5),
linewidth = 0.95,
transform = ccrs.PlateCarree(),
zorder = 20,
regrid_shape = 15,
color = 'tab:blue')
elif draw == "quiver": # 矢印
cs = ax.contour(lons, lats, data[::-1, :], levels = levels, linewidths = 3)
ax.clabel(cs, fmt = '%d', fontsize = 14)
ax.quiver(lons, lats, u, v, transform = ccrs.PlateCarree(), regrid_shape = 20, color='tab:red')
plt.show()
In [6]:
# gzipされたファイルを解凍し、バイナリデータを読み込む
def Read_gz_Binary(file):
file_tmp = file + "_tmp"
with gzip.open(file, 'rb') as f_in:
with open(file_tmp, 'wb') as f_out:
shutil.copyfileobj(f_in, f_out)
bin_data = np.fromfile(file_tmp, np.float32)
os.remove(file_tmp)
return bin_data.reshape( [168, 128] )
In [7]:
dir2 = "./signate/met_data_2017_01/"
# 850 hpa気温（850 hPaの高度は一定値ではないがおよそ1500m)
file0411 = dir2 + "train/met/2017/04/11/TMP.850.3.2017041100.gz"
data0411 = Read_gz_Binary( file0411 )
draw_weather_map(data0411-273.15, draw = "shaded", levels = list(np.arange(-15, 25, 3))) # 温度３度おき
In [8]:
# 欠損データ穴埋め処理
def fill_lack_data(data):
# まず欠損行の穴埋めは、値が存在する上下端の行の値をそのままコピーする
data[0:2] = data[2]
data[154:] = data[153]
# 欠損列の穴埋めも、値が存在する左右端の列の値をそのままコピーする
data[:, :8] = data[:, 8].reshape(-1, 1)
return data
In [9]:
# 欠損箇所の簡易埋め
dir2 = "./signate/met_data_2017_01/"
file0411 = dir2 + "train/met/2017/04/11/TMP.850.3.2017041100.gz"
data0411 = Read_gz_Binary( file0411 )
data0411 = fill_lack_data( data0411 )
draw_weather_map(data0411-273.15, draw = "shaded", levels = list(np.arange(-15, 25, 3))) # 温度３度おき
In [10]:
# 海面気圧 等圧線4hPa
dir2 = "./signate/met_data_2017_01/"
file0411 = dir2 + "train/met/2017/04/11/PRMSL.msl.3.2017041100.gz"
data0411 = Read_gz_Binary( file0411 )
data0411 = fill_lack_data( data0411 )
draw_weather_map(data0411/100, draw = "contour", levels = list(np.arange(992, 1028, 4))) # 等圧線4hPa
In [11]:
# 700 hPaの湿度
dir2 = "./signate/met_data_2017_01/"
file0411 = dir2 + "train/met/2017/04/11/RH.700.3.2017041100.gz"
data0411 = Read_gz_Binary( file0411 )
data0411 = fill_lack_data( data0411 )
draw_weather_map(data0411, draw = "shaded", levels = list(np.arange(0, 101, 5))) # 5%間隔
In [12]:
# 500 hPaの高度
dir2 = "./signate/met_data_2017_01/"
file0411 = dir2 + "train/met/2017/04/11/HGT.500.3.2017041100.gz"
data0411 = Read_gz_Binary( file0411 )
data0411 = fill_lack_data( data0411 )
draw_weather_map(data0411, draw = "contour", levels = list(np.arange(5220, 5880, 60))) # 60m間隔
In [13]:
# 300 hPaの風
dir2 = "./signate/met_data_2017_01/"
file0411a = dir2 + "train/met/2017/04/11/UGRD.300.3.2017041100.gz" # u wind 東西風 4月11日
file0411b = dir2 + "train/met/2017/04/11/VGRD.300.3.2017041100.gz" # w wind 南北風
file0411c = dir2 + "train/met/2017/04/11/HGT.300.3.2017041100.gz" # 等高線
data0411a = Read_gz_Binary( file0411a )
data0411a = fill_lack_data( data0411a )
data0411b = Read_gz_Binary( file0411b )
data0411b = fill_lack_data( data0411b )
data0411c = Read_gz_Binary( file0411c )
data0411c = fill_lack_data( data0411c )
draw_weather_map(data0411c,
u = data0411a/0.5144,
v = data0411b/0.5144,
draw = "barbs",
levels = list(np.arange(8400, 9840, 120))) #
In [14]:
draw_weather_map(data0411c,
u = data0411a/0.5144,
v = data0411b/0.5144,
draw = "quiver",
levels = list(np.arange(8400, 9840, 120))) #
In [15]:
# 風の東西成分(u)・南北成分(v)から、風向風速を計算する
def Calc_Wind(u, v):
wspd = np.sqrt( u**2 + v**2 )
wdir = np.rad2deg( np.arctan2(u, v) )+180
return np.round(wspd, 1), np.round(wdir, 0)
In [16]:
Calc_Wind(u = 2, v = -3) # 西風2m/s 北風3m/s 結果 風速、風向
Out[16]:
In [17]:
#評価領域の切り取り
valid_img = img4[40:40+420, 130:130+340]
print( "image shape : ", valid_img.shape )
print( "data type : ", valid_img.dtype )
fig, axes = plt.subplots(1, 1, figsize = (12, 8))
axes.imshow(valid_img, cmap = 'gray')
axes.set_title("SatelliteImage on 2017/04/01 00:00UT")
plt.show()
In [ ]:
In [ ]: