[PYTHON] Introduction of drawing code for figures with a certain degree of perfection of meteorological data
This article is the 11th day of Student LT Advent Calendar 2019 --Qiita.
Hello. I am a graduate student majoring in atmospheric science.
This time, as you can see in the previous article, I focused on drawing weather data using Python. The drawing code to be introduced is available on GitHub, so if you are interested, please download it.
Introduction
The following articles were found in relation to the drawing and handling of meteorological data.
It is very easy to understand, but both are just an introduction to drawing one element.
I would like to differentiate from the above article by introducing the code from organizing the data of the figure with a certain degree of perfection to drawing as shown below.
Click here for the rendering to be created this time. The target captures the case of Typhoon Lionrock 2016.
| ・ Fig. 1 | ・ Fig. 2 |
|---|---|
 |
 |
Introduction of each figure
I will introduce my own interpretation of the figure. (... demand may be low)
- Figure 1
We are drawing a combination of meteorological elements from 06:00 on August 21, 2016 to 00:00 on August 31, 2016. The contour represents the altitude field, the blue vector represents the wind speed field of 10 m / s or more, and the diagonal line + yellow part represents the particularly strong cyclone part in the specified altitude field.
- I decided the extraction method of the strong low pressure part at my own discretion this time, but I am not deeply elaborate because it is not academic but an introduction of drawing method.
→ You can see how the power is gradually expanding before landing.
- Figure 2
Hexagrid expresses which point the typhoon has passed through for the last 20 years (2000-2019). In addition, the plot shows a specific typhoon path, and this time the path of typhoon No. 10 in 2016 is plotted according to Fig. 1.
→ It can be read from the hexagrid that most of the typhoons so far have passed near the eastern side of Taiwan and moved northward and approached Japan, but it was confirmed that typhoon No. 10 in 2016 was a special route. I can do it.
Data to prepare
- Global Numerical Weather Prediction Model GPV
It can be obtained from Research Institute for Sustainable Humanosphere (RISH). I will. Please see Weather Line Support Center for detailed grid information and delivery time.
- Please note the following points when using the above site.
We provide data for educational and research institutes. If you frequently need data for corporate activities, please purchase the data directly from the Meteorological Business Support Center and cooperate in the maintenance and development of the entire data provision scheme. (Quoted from the site)
The data used this time is the data obtained from the above site, and the data of east-west wind, north-south wind, vertical wind, temperature, and altitude field are cut out in binary format using wgirb2.
The cut out data format is described in Fortran90. ... for reference only.
format.f90
integer, parameter :: nx = 720, ny = 361, np = 12
integer, parameter :: elem = 5
real(4) :: gpv(elem, hgt, nx, ny)
real(4) :: pr(hgt)
character(4) :: elem_list(elem)
data pr / 1000., 925., 850., 700., 600., 500., 400., 300., 250., 200., 150., 100. /
data elem_list / "U ", "V ", "W ", "T ", "Z " /
open(*, fille=*, form='unformatted', access='direct', recl=4*nx*ny*np*elem)
close(*)
- Typhoon location table (csv format)
It can be obtained by year from Japan Meteorological Agency. Information from the occurrence of a typhoon to an extratropical cyclone is provided at 6-hour intervals.
Execution code
I hope you read both codes based on main_driver. If you want to execute it, please change the PATH setting.
Meteorological element drawing code (Fig. 1)
First, I will introduce the code that draws the weather element.
weather_bin2plot.py
# -*- coding: utf-8 -*-
"""
Created on 2019.12.11
@author: Toyo_Daichi
"""
import numpy as np
import pandas as pd
import glob
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
class mapping:
def __init__(self):
pass
def gpv_data_coef(self):
nx, ny, hgt = 720, 361, 12
elem = 5
return nx, ny, hgt, elem
def prj_coef(self, nx, ny):
dx, dy = 0.5, 0.5
lon, lat = [], []
for ix in range(nx):
lon += [ float('{:.2f}'.format(dx*ix)) ]
for iy in range(ny):
lat += [ float('{:.2f}'.format(-90. + dy*iy)) ]
X, Y = np.meshgrid(lon, lat)
return X, Y
def preparating_data(self, yyyy, mm, dd, hh):
time_list = []
num_time_list = len(time_list)
month_thirtyone = [ 1, 3, 5, 7, 8, 10, 12 ]
month_thirty = [ 4, 6, 9, 11 ]
month_twntynine = [ 2 ]
while num_time_list < 40:
time_list.append(str(yyyy) + str('%02d' % mm) + str('%02d' % dd) + str('%02d' % hh) + '00')
hh = hh - 6
if hh < 0 and dd == 1:
mm, hh = mm - 1, 18
if mm in month_thirty:
dd = 30
elif mm in month_thirtyone:
dd = 31
elif mm in month_twntynine:
if yyyy % 4 == 0:
dd = 28
else:
dd =29
elif hh < 0:
dd, hh = dd - 1, 18
num_time_list += 1
return time_list
def setup_gpv_filelist(self, path, time_list):
filelist = []
for i_file in time_list:
filelist.append(path + '/data')
return filelist
def open_gpv_filelist(self, gpv_filelist, nx, ny, hgt, elem):
data = [ []*i for i in range(len(gpv_filelist)) ]
for num_gpv, gpv_file in enumerate(gpv_filelist):
with open(gpv_file, 'rb') as ifile:
"""
elem: 1-u, 2-v, 3-w, 4-tmp, 5-height
"""
data[num_gpv] = np.fromfile(ifile, dtype='>f', sep = '').reshape(elem, hgt, ny, nx)
print('..... Preparating data for ' + gpv_file, ', shapes :: ', data[num_gpv].shape)
return data
def main_mapping_tool(self, mode, path, time_list, nx, ny, *, gpv_datalist='None', snap_step=0, level='1000'):
fig, ax = plt.subplots()
outpath = path + '/fig/'
lat_min, lat_max = 17, 50
lon_min, lon_max = 120, 155
mapping = Basemap( projection='lcc', resolution="i", lat_0=35, lon_0=140, fix_aspect=(1,1), llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max )
mapping.drawcoastlines(color='black', linewidth=0.5)
mapping.drawmeridians(np.arange(0, 360, 5), labels=[False, True, False, True], fontsize='small', color='gray', linewidth=0.5)
mapping.drawparallels(np.arange(-90, 90, 5), labels=[True, False, False, True], fontsize='small', color='gray', linewidth=0.5)
lon_list, lat_list = self.prj_coef(nx, ny)
x, y = mapping(lon_list, lat_list)
if mode == 1 or mode == 2:
if level == "1000":
hPa, min_list = 0, [0, 200]
elif level == "925":
hPa, min_list = 1, [0, 600]
elif level == "850":
hPa, min_list = 2, [0, 750]
elif level == "700":
hPa, min_list = 3, [0, 1500]
elif level == "600":
hPa, min_list = 4, "None"
elif level == "500":
hPa, min_list = 5, "None"
elif level == "400":
hPa, min_list = 6, "None"
elif level == "300":
hPa, min_list = 7, "None"
elif level == "250":
hPa, min_list = 8, "None"
elif level == "200":
hPa, min_list = 9, "None"
elif level == "150":
hPa, min_list = 10, "None"
elif level == "100":
hPa, min_list = 11, "None"
gpv_u_data = gpv_datalist[snap_step][0][hPa]
gpv_v_data = gpv_datalist[snap_step][1][hPa]
speed = np.sqrt(gpv_u_data*gpv_u_data + gpv_v_data*gpv_v_data)
speed_list = list(range(0, 50, 25))
gpv_height_data = gpv_datalist[snap_step][4][hPa]
num_list = list(range(0, 7500, 10))
contour = mapping.contour(x, y, gpv_height_data, linewidths=0.25, linestyles='-', levels=num_list, colors='m')
contour.clabel(fmt='%1.1f', fontsize=6.5)
if not min_list == "None":
lines = mapping.contourf(x, y, gpv_height_data, min_list, alpha=0.5, hatches=['///'], lw=1., zorder=0, edgecolor='r', colors="y")
for i_nx, i_ny in itertools.product(range(nx), range(ny)):
if speed[i_ny][i_nx] > 10 and lon_min <= lon_list[i_ny][i_nx] <= lon_max and lat_min <= lat_list[i_ny][i_nx] <= lat_max:
print('...... Halfway step, ', i_nx, i_ny, speed[i_ny][i_nx])
vector = mapping.quiver(x[i_ny][i_nx], y[i_ny][i_nx], gpv_u_data[i_ny][i_nx], gpv_v_data[i_ny][i_nx], color='c', units='dots', scale=2.0, alpha=0.6)
plt.title(time_list[snap_step] + ' @GSM ' + level + 'hPa' , loc='left', fontsize=10)
plt.quiverkey(vector, 0.75, 0.9, 10, '10 m/s', labelpos='W', coordinates='figure')
plt.savefig(outpath + 'GPV_elem_' + time_list[snap_step] + '.png')
print('...... Saving fig :: ', outpath + 'GPV_elem_' + time_list[snap_step] + '.png')
#plt.show()
def main_driver(self, mode, indir, time_list, level):
nx, ny, hgt, elem = self.gpv_data_coef()
gpv_filelist = self.setup_gpv_filelist(indir, time_list)
gpv_datalist = self.open_gpv_filelist(gpv_filelist, nx, ny, hgt, elem)
if mode == 1:
snap_step = 0
self.main_mapping_tool(mode, indir, time_list, nx, ny, gpv_datalist=gpv_datalist, level=level, snap_step=snap_step)
elif mode == 2:
for snap_step in range(len(gpv_datalist)):
self.main_mapping_tool(mode, indir, time_list, nx, ny, gpv_datalist=gpv_datalist, level=level, snap_step=snap_step)
if __name__ == "__main__":
mapp = mapping()
#input dir
indir = '/your_directory/'
#target time
yyyy, mm, dd, hh = 2016, 8, 31, 0
time_list = mapp.preparating_data(yyyy, mm, dd, hh)
"""
Target altitude list
1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100 hPa
"""
level = '925'
#choose mode, if you append new func. more anynum.
"""
2019.12.11
mode 1: Normal weather element info at GPV GSM snap shot.
mode 2: Normal weather element info at GPV GSM gif.
"""
mode = 1
#main_driver
mapp.main_driver(mode, indir, time_list, level)
The points I devised with this code were as follows.
- Draw the wind speed vector only at 10m / s or more.
- Draw large values at each altitude with diagonal lines.
- Allows you to select between gif format image creation and snapshot image creation in mode.
Typhoon track drawing code (Fig. 2)
Next, I will introduce the degree of typhoon course for the past 20 years and the code that created the course of a specific typhoon.
typhoon_csv2plot.py
# -*- coding: utf-8 -*-
"""
Created on 2019.12.11
@author: Toyo_Daichi
"""
import numpy as np
import pandas as pd
import glob
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
class mapping:
def __init__(self):
pass
def setup_csv_filelist(self, path, *, year='*'):
filelist = []
filelist = glob.glob(path + 'table' + str(year) + '.csv')
return filelist
def open_csv_filelist(self, csv_filelist, *, typhoon_number='None'):
csv_datalist = [ []*i for i in range(len(csv_filelist)) ]
list_num = list(range(11))
list_option = ( 'year', 'month', 'day', 'hour(UTC)', 'typhoon_number', 'typhoon_name', 'rank','latitude', 'longitude', 'central_pressure', 'max_wind')
for num_file, infile in enumerate(csv_filelist):
print('..... Preparating data for ' + str(num_file) + ' ' + str(infile))
tmp_data = pd.read_csv(infile, usecols=list_num, skiprows=1, names=list_option, sep=',')
if typhoon_number == "None":
tmp_lat_list = tmp_data['latitude'].values.tolist()
csv_datalist[num_file].append(tmp_lat_list)
tmp_lon_list = tmp_data['longitude'].values.tolist()
csv_datalist[num_file].append(tmp_lon_list)
tmp_centpre_list = tmp_data['central_pressure'].values.tolist()
csv_datalist[num_file].append(tmp_centpre_list)
else:
csv_datalist = []
specific_data = tmp_data.query('typhoon_number == %s' % typhoon_number)
tmp_lat_list = specific_data['latitude'].values.tolist()
csv_datalist.append(tmp_lat_list)
tmp_lon_list = specific_data['longitude'].values.tolist()
csv_datalist.append(tmp_lon_list)
tmp_centpre_list = specific_data['central_pressure'].values.tolist()
csv_datalist.append(tmp_centpre_list)
return csv_datalist
def main_mapping_tool(self, path, csv_datalist, csv_specific_datalist='None', typhoon_info='None'):
fig, ax = plt.subplots()
outpath = path + '/fig/'
lat_min, lat_max = 17, 50
lon_min, lon_max = 120, 155
mapping = Basemap( projection='lcc', resolution="i", lat_0=35, lon_0=140, fix_aspect=(1,1), llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max )
mapping.drawcoastlines(color='black', linewidth=0.5)
mapping.drawmeridians(np.arange(0, 360, 5), labels=[False, True, False, True], fontsize='small', color='gray', linewidth=0.5)
mapping.drawparallels(np.arange(-90, 90, 5), labels=[True, False, False, True], fontsize='small', color='gray', linewidth=0.5)
full_lat_list = list(map(lambda x: x[0], csv_datalist))
full_lon_list = list(map(lambda x: x[1], csv_datalist))
full_lat_list = np.sum(full_lat_list, axis=0)
full_lon_list = np.sum(full_lon_list, axis=0)
lat_list, lon_list = [], []
for i_num in range(len(full_lat_list)):
if(lat_min <= full_lat_list[i_num] <= lat_max and lon_min <= full_lon_list[i_num] <= lon_max):
lat_list.append(full_lat_list[i_num])
lon_list.append(full_lon_list[i_num])
x, y = mapping(lon_list, lat_list)
hexbin = mapping.hexbin(np.array(x), np.array(y), gridsize=[30, 30], cmap='Blues', edgecolors='gray', mincnt=8)
if not csv_specific_datalist == "None":
specific_lat_list, specific_lon_list = csv_specific_datalist[0], csv_specific_datalist[1]
specific_centpre_list = csv_specific_datalist[2]
lat_list, lon_list, centpre_list = [], [], []
for i_num in range(len(specific_lon_list)):
if(lat_min <= specific_lat_list[i_num] <= lat_max and lon_min <= specific_lon_list[i_num] <= lon_max):
lat_list.append(specific_lat_list[i_num])
lon_list.append(specific_lon_list[i_num])
centpre_list.append(specific_centpre_list[i_num])
case_x, case_y = mapping(lon_list, lat_list)
mapping.plot(case_x, case_y, linewidth=0.5, color='red', ls='--', marker='o', ms=2)
cbar = plt.colorbar(hexbin)
cbar.set_label('number', fontsize=8)
if not csv_specific_datalist == "None":
plt.title('Course of typhoon 2000-2019' + ', ' + 'Typhoon track: ' + str(typhoon_info[1]), loc='left', fontsize=9)
else:
plt.title('Course of typhoon 2000-2019', loc='left', fontsize=10)
plt.savefig(outpath + 'typhoon_strength2000-2019' + '_track_' + str(typhoon_info[1]) + '.png')
print('..... savefig ' + outpath + 'typhoon_strength2000-2019' + '_track_' + str(typhoon_info[1]) + '.png')
#plt.show()
def main_driver(self, indir, *, typhoon_info='None'):
csv_filelist = self.setup_csv_filelist(indir)
csv_datalist = self.open_csv_filelist(csv_filelist)
"""
csv_datalist.shape = [year][latitude][longitude][central_pressure]
"""
if not typhoon_info == "None":
print('..... Check specific case filelist')
csv_specific_filelist = self.setup_csv_filelist(indir, year=typhoon_info[0])
csv_specific_datalist = self.open_csv_filelist(csv_specific_filelist, typhoon_number=typhoon_info[1])
self.main_mapping_tool(indir, csv_datalist, csv_specific_datalist=csv_specific_datalist, typhoon_info=typhoon_info)
else:
self.main_mapping_tool(indir, csv_datalist)
if __name__ == "__main__":
mapp = mapping()
#input dir
indir = '/your_directory/'
"""
If you want to write a specific typhoon route, enter the typhoon information.
typhoon_info = [year, typhoon_number]
"""
typhoon_info = [2016, 1610]
#main_driver
mapp.main_driver(indir, typhoon_info=typhoon_info)
print('Normal END')
The points I devised with this code were as follows.
- Draw the frequency of typhoon paths for 20 years with a hexagrid.
- By describing a specific typhoon number, the route of that typhoon is drawn.
Usage environment used this time
I'm sorry to describe the environment.
- Python 3.6.7
requirement.txt
basemap==1.2.0
matplotlib==3.1.1
Finally
The code I created is small, but I've made some ingenuity in my own way, so I'd like to introduce it in the future. If you find out how to write better, please let us know in the comments.
Thank you for reading to the end.
Referenced site
- Meteorology x Python ~ Satellite precipitation data ~ https://qiita.com/OSAKO/items/d54da432fff34d4be1d8
Recommended Posts