[PYTHON] Einführung des Zeichnungscodes für Figuren mit einem gewissen Grad an Perfektion der Wetterdaten
Dieser Artikel ist der 11. Tag von Student LT Adventskalender 2019 - Qiita.
Hallo. Ich bin ein Doktorand mit Schwerpunkt Atmosphärenwissenschaften.
Wie Sie im vorherigen Artikel sehen können, habe ich mich dieses Mal darauf konzentriert, Wetterdaten mit Python zu zeichnen. Der einzuführende Zeichnungscode ist auf GitHub verfügbar. Wenn Sie also interessiert sind, laden Sie ihn bitte herunter.
Einführung
Die folgenden Artikel beziehen sich auf das Zeichnen und Behandeln von meteorologischen Daten.
- Visualisieren Sie GRIB2 auf einer Karte mit Cartopy
- Visualisiere grib2 auf einer Karte mit Python (matplotlib)
Es ist sehr leicht zu verstehen, aber beide sind nur eine Einführung in das Zeichnen eines Elements.
Ich möchte mich von dem obigen Artikel unterscheiden, indem ich den Code von der Organisation der Daten der Figur mit einem gewissen Grad an Perfektion bis zur Zeichnung wie unten gezeigt einführe.
Klicken Sie hier, um die fertige Zeichnung anzuzeigen, die dieses Mal erstellt werden soll. Das Ziel erfasst den Fall des Taifuns Nr. 10 im Jahr 2016.
| ・ Abb. 1 | ・ Abb. 2 |
|---|---|
 |
 |
Einführung jeder Figur
Ich werde meine eigene Interpretation der Figur vorstellen. (... Nachfrage kann gering sein)
- Abbildung 1
Wir zeichnen eine Kombination meteorologischer Elemente von 06:00 Uhr am 21. August 2016 bis 00:00 Uhr am 31. August 2016. Die Kontur repräsentiert das Höhenfeld, der blaue Vektor repräsentiert das Windgeschwindigkeitsfeld von 10 m / s oder mehr und die diagonale Linie + der gelbe Teil repräsentiert den besonders starken Niederdruckteil im angegebenen Höhenfeld.
- Ich habe mich diesmal nach eigenem Ermessen für die Methode zum Extrahieren des starken Niederdruckteils entschieden, bin aber nicht sehr ausführlich, da es sich nicht um eine akademische Methode handelt, sondern um eine Einführung in die Zeichenmethode.
→ Sie können sehen, wie sich die Leistung vor der Landung allmählich erweitert.
- Figur 2
Hexagrid drückt aus, welche Punkte die Taifune der letzten 20 Jahre (2000-2019) durchlaufen haben. Darüber hinaus zeigt das Diagramm einen bestimmten Taifunpfad, und diesmal ist der Pfad des Taifuns Nr. 10 im Jahr 2016 gemäß Abb. 1 dargestellt.
→ Aus dem Hexagrid geht hervor, dass die meisten Taifune bisher in der Nähe der Ostseite Taiwans vorbeigefahren sind und sich nach Norden bewegt haben und sich Japan näherten. Es wurde jedoch bestätigt, dass der Taifun Nr. 10 im Jahr 2016 eine Sonderroute war. Ich kann es schaffen
Daten vorzubereiten
- Globales numerisches Prognosemodell GPV
Es kann von Forschungsinstitut für nachhaltige Humanosphäre (RISH) bezogen werden. Ich werde. Ausführliche Informationen zum Raster und zur Lieferzeit finden Sie im Meteorological Line Support Center.
- Bitte beachten Sie die folgenden Punkte, wenn Sie die oben genannte Website verwenden.
Wir liefern Daten für Bildungs- und Forschungseinrichtungen. Wenn Sie häufig Daten für Unternehmensaktivitäten benötigen, kaufen Sie diese bitte direkt beim Meteorological Business Support Center und kooperieren Sie bei der Wartung und Entwicklung des gesamten Datenbereitstellungsschemas. (Zitiert von der Website)
Die Daten, die dieses Mal verwendet werden, sind die Daten, die von der obigen Stelle erhalten wurden, und die Daten von Ost-West-Wind, Nord-Süd-Wind, vertikalem Wind, Temperatur und Höhenfeld werden unter Verwendung von wgirb2 im Binärformat ausgeschnitten.
Das ausgeschnittene Datenformat ist in Fortran90 beschrieben. ... nur als Referenz.
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(*)
- Taifun-Standorttabelle (CSV-Format)
Es kann nach Jahr bei Meteorological Agency bezogen werden. Informationen vom Zeitpunkt des Taifuns bis zum gemäßigten Niederdruck werden in Intervallen von 6 Stunden beschrieben.
Ausführungscode
Ich hoffe, Sie lesen beide Codes basierend auf main_driver. Wenn Sie es ausführen möchten, ändern Sie bitte die PATH-Einstellung.
Zeichnungscode für meteorologische Elemente (Abb. 1)
Zuerst werde ich den Code einführen, der die meteorologischen Elemente zeichnet.
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)
Die von diesem Code entwickelten Punkte waren wie folgt.
- Zeichnen Sie den Windgeschwindigkeitsvektor nur mit 10 m / s oder mehr.
- Zeichnen Sie in jeder Höhe große Werte mit diagonalen Linien.
- Ermöglicht die Auswahl zwischen der Erstellung von Bildern im GIF-Format und der Erstellung von Schnappschussbildern im Modus.
Zeichnungscode für die Taifunspur (Abb. 2)
Als nächstes werde ich den Grad des Taifunkurses der letzten 20 Jahre und den Code vorstellen, der den Verlauf eines bestimmten Taifuns erstellt hat.
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')
Die von diesem Code entwickelten Punkte waren wie folgt.
- Zeichnen Sie die Häufigkeit von Taifunpfaden für 20 Jahre mit einem Hexagrid.
- Durch die Beschreibung einer bestimmten Taifunnummer wird die Route dieses Taifuns gezeichnet.
Diesmal verwendete Nutzungsumgebung
Es tut mir leid, die Umgebung zu beschreiben.
- Python 3.6.7
requirement.txt
basemap==1.2.0
matplotlib==3.1.1
Schließlich
Der Code, den ich erstellt habe, ist klein, aber ich habe auf meine eigene Weise etwas Einfallsreichtum entwickelt, sodass ich ihn in Zukunft gerne einführen möchte. Wenn Sie herausfinden, wie Sie besser schreiben können, würde ich mich freuen, wenn Sie mich in den Kommentaren unterrichten könnten.
Vielen Dank für das Lesen bis zum Ende.
Referenzierte Site
- Meteorologie x Python ~ Satellitenniederschlagsdaten ~ https://qiita.com/OSAKO/items/d54da432fff34d4be1d8
Recommended Posts