[PYTHON] Introduction du code de dessin pour les figures avec un certain degré de perfection des données météorologiques
Cet article est le 11ème jour du Calendrier de l'Avent étudiant LT 2019 - Qiita.
Bonjour. Je suis un étudiant diplômé spécialisé en sciences atmosphériques.
Cette fois, comme vous pouvez le voir dans l'article précédent, je me suis concentré sur le dessin de données météorologiques à l'aide de Python. Le code de dessin à introduire est disponible sur GitHub, donc si vous êtes intéressé, veuillez le télécharger.
introduction
Les articles suivants ont été trouvés liés au dessin et à la gestion des données météorologiques.
C'est très facile à comprendre, mais les deux ne sont qu'une introduction au dessin d'un élément.
Je voudrais me différencier de l'article ci-dessus en introduisant le code de l'organisation des données de la figure avec un certain degré de perfection au dessin comme indiqué ci-dessous.
Cliquez ici pour que le dessin terminé soit créé cette fois. La cible capture le cas du typhon n ° 10 en 2016.
| · Fig. 1 | ・ Fig. 2 |
|---|---|
 |
 |
Présentation de chaque figure
Je vais vous présenter ma propre interprétation de la figure. (... la demande peut être faible)
- Figure 1
Nous dessinons une combinaison d'éléments météorologiques de 06h00 le 21 août 2016 à 00h00 le 31 août 2016. Le contour représente le champ d'altitude, le vecteur bleu représente le champ de vitesse du vent de 10 m / s ou plus, et la ligne diagonale + partie jaune représente la partie basse pression particulièrement forte dans le champ d'altitude spécifié.
- J'ai décidé cette fois la méthode d'extraction de la partie forte basse pression à ma discrétion, mais je ne suis pas profondément élaborée car ce n'est pas académique mais une introduction à la méthode de dessin.
→ Vous pouvez voir comment la puissance augmente progressivement avant d'atterrir.
- Figure 2
Hexagrid exprime quels points les typhons des 20 dernières années (2000-2019) sont passés. De plus, le graphique montre une trajectoire de typhon spécifique, et cette fois la trajectoire du typhon n ° 10 en 2016 est tracée selon la Fig.1.
→ On peut lire à partir de l'hexagride que la plupart des typhons jusqu'à présent sont passés près du côté est de Taiwan et se sont déplacés vers le nord et ont approché le Japon, mais il a été confirmé que le typhon n ° 10 en 2016 était une route spéciale. Je peux le faire.
Données à préparer
- Modèle de prévision numérique global GPV
Il peut être obtenu auprès de Institut de recherche pour une humanosphère durable (RISH). Je vais. Veuillez consulter le Centre d'assistance de la ligne météorologique pour des informations détaillées sur la grille et le délai de livraison.
- Veuillez noter les points suivants lors de l'utilisation du site ci-dessus.
Nous fournissons des données aux établissements d'enseignement et de recherche. Si vous avez fréquemment besoin de données pour les activités de l'entreprise, veuillez acheter les données directement auprès du Centre d'assistance météorologique aux entreprises et coopérer à la maintenance et au développement de l'ensemble du système de fourniture de données. (Cité sur le site)
Les données utilisées cette fois sont les données obtenues à partir du site ci-dessus, et les données de vent est-ouest, vent nord-sud, vent vertical, température et champ d'altitude sont découpées au format binaire à l'aide de wgirb2.
Le format des données découpées est décrit dans Fortran90. ... pour référence seulement.
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(*)
- Table de localisation Typhoon (format csv)
Il peut être obtenu par année auprès de Agence météorologique. Les informations depuis le moment du typhon jusqu'à la dépression tempérée sont décrites toutes les 6 heures.
Code d'exécution
J'espère que vous avez lu les deux codes basés sur main_driver. Si vous souhaitez l'exécuter, veuillez modifier le paramètre PATH.
Code de dessin des éléments météorologiques (Fig.1)
Tout d'abord, je présenterai le code qui dessine les éléments météorologiques.
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)
Les points élaborés par ce code étaient les suivants.
- Dessinez le vecteur de vitesse du vent uniquement à 10 m / s ou plus.
- Dessinez de grandes valeurs à chaque altitude avec des lignes diagonales.
- Vous permet de choisir entre la création d'image au format gif et la création d'image instantanée en mode.
Code de dessin de la trajectoire du typhon (Fig.2)
Ensuite, je présenterai le degré de cours de typhon au cours des 20 dernières années et le code qui a créé le cours d'un typhon spécifique.
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')
Les points élaborés par ce code étaient les suivants.
- Dessinez la fréquence des trajectoires des typhons pendant 20 ans avec un hexagrid.
- En décrivant un numéro de typhon spécifique, l'itinéraire de ce typhon est tracé.
Environnement d'utilisation utilisé cette fois
Je suis désolé de décrire l'environnement.
- Python 3.6.7
requirement.txt
basemap==1.2.0
matplotlib==3.1.1
finalement
Le code que j'ai créé est petit, mais j'ai fait preuve d'ingéniosité à ma manière, je voudrais donc le présenter à l'avenir. Si vous apprenez à mieux écrire, je vous serais reconnaissant de bien vouloir m'apprendre dans les commentaires.
Merci d'avoir lu jusqu'au bout.
Site référencé
- Météorologie x Python ~ Données de précipitations satellites ~ https://qiita.com/OSAKO/items/d54da432fff34d4be1d8
Recommended Posts